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Abstract In this lectures various methods which give a possibility 
to extend an area of applicability of perturbation series and hence to 
omit their local character are analysed. While applying asymptotic 
methods as a rule the following situation appears: the existence of 
asymptotics for e — >• implies an existence of the asymptotics for 
e — >• oo. Therefore, the idea to construct one function valid for the 
whole parameter interval for e is very attractive. The construction 
of asymptotically equivalent functions possessing a known asymp- 
totic behaviour for e — >• and e — >• oo will be discussed. Using 
summation and interpolation procedures we focus on continuous 
models derived from a discrete micro- structure. Various contin- 
ualization procedures that take the non-local interaction between 
variables of the discrete media into account are analysed. 
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1 Overview of Manuscript Preparation and Delivery 

Asymptotic analysis is a constantly growing branch of mathematics which 
influences the de velopment of vari ous pure an d applied sciences. The famous 
mathematicians iFriedrichsl ([1955) and lSegell (l966) said that an asymptotic 



description is not only a suitable instrument for the mathematical analysis 
of nature but that it also has an additional deeper intrinsic meaning, and 
that the asymptotic approach is more than just a mathematical technique; 
it plays a rather fundamental role in science. And here it appears that 
the many existing asymptotic methods comprise a set of appro aches that 
in some way belong rather to art than to science. iKruskall ([19631 ) even in- 



troduced special term "asymptotology" and defined it as the "art of dealing 
with applied mathematical systems in limiting cases". Here it should be 
noted that he called for a formalization of the accumulated experience to 
convert the art of asymptotology into a science of asymptotology. 
Asymptotic methods for solving mechanical and physical problems have 
been deve loped by man y authors. We can mention e d the excellent mono- 
graphs bv [Hind] (Il99lh ["Kevorkian andColel (fl996h . iMillerl (f2QQ6h . iNavfehl 



(jl98ll . l2QQ0h . Ivan Dvkd jl975a.b). VerhuTstl 



2005), and many others. The 



main feature of the present book can be formulated as follows: it deals 
with new trends and applications of asymptotic approaches in the fields 
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of Nonlinear Mechanics and Mechanics of Solids. It illuminates the devel- 
opments in the field of asymptotic mathematics from different viewpoints, 
reflecting the field's multidisciplinary nature. The choice of topics reflects 
the authors' own research experience and participation in applications. The 
authors have paid special attention to examples and discussions of results, 
and have tried to avoid burying the central ideas in formalism, notations, 
and technical details. 



2 Some Nonstandard Perturbation Procedures 

2.1 Choice of Small Parameters 

The choice of the asymptotic method and the introduction of small di- 
mensionless parameters in the system is very often the most significant and 
informal part of the analytical study of physical problems. This should 
help the experience and intuition, analysis of the physical nature of the 
problem, experimental and numerical results. It often dictated by physical 
considerations, is clearly shown as dimensionless and scaling. However, it is 
sometimes advantageous to use is not obvious, and perhaps even strange at 
first glan ce, the initial approx imation. To illustrate this, consider a simple 



example ([Bender et al.L 119981 ) : algebraic equation 

x 5 +x = l. (1) 



We seek the real root of Eq. (pQ) , the exact value of which can be determined 
numerically: x — 0.75487767... . A small parameter e is not included 
explicitly in Eq. (pQ). Consider the various possibilities of introducing a 
parameter e in Eq. (pQ). 

1. We introduce a small parameter e with a nonlinear term of Eq. (pQ) 

ex 5 +x = l, (2) 



and present x as a series of e 

x = a + a\s + a 2 e 2 + (3) 



Substituting the series (|3]) into Eq. (|2]) and equating terms of equal 
powers, we obtain 

ciq = 1, cli = — 1, (i2 = 5, a 3 = —35, a4 = 285, a 5 = —2530, 
a 6 = 23751. 
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These values can be obtained by a closed expression for the coeffi- 
cients a n : 

(-l) n (5n)! 
an ~ n!(4n + l)f 

The radius R of convergence of the series (|3]) is R = |s = 0.08192. 
Consequently, for e = 1 series (|3]) diverges very fast, so the sum of 
the first six terms is 21476. The situation can be corrected by the 
method of Pade approximants (see Sect. 2). Constructing a Pade 
approximant with three terms in the numerator and denominator and 
calculating it with e = 1, we obtain the value of the root x = 0.76369 
(the deviation from the exact value is 1.2%). 
2. We now introduce a small parameter e near the linear term of Eq. (pQ) 

£ 5 +££=l. (4) 

Presenting the solution of Eq. (HJ in the form 

x(e) = b + he + b 2 e 2 + . . . , (5) 
we have, after applying the standard procedure of perturbation method 



b = 1, h = -1, b 2 = -I25, «3 = -I125, h = 0, 65 = 

21 h = 78 Zb iZb 

15625' ^ 6 78125* 

And in this case we can construct a general expression for the co- 
efficients with: 

r [(An - l)/5] 
5r[(4-n)/5]n!' 

and determine the radius of convergence of the series R = 
4 (4/ 5 ) = 1, 64938 .... The value of x(l), taking into account the first 
six terms of the series (j5j), deviates from the exact by 0.07%. 
3. We introduce the "small parameter" 5 in the exponent 

x 1+5 +£ = l, (6) 

and represent x in the form 

x = c + c\ 5 + c 2 £ 2 + (7) 

In addition, we use the expansion: 

= x(l + *ln|x| + ...). 
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Coefficients of the series ([7]) are determined easily: 



c =0.5, ci=0.251n2, c 2 



-0.125 In 2, 



The radius of convergence is one in this case, and calculate when it 
should be S = 4. Using Pade approximants with three terms in the nu- 
merator and denominator, if e — 1, we find x = 0.75448 that only devi- 
ates from the exact by 0.05% result. Calculating q for i = 0, 1, . . . , 12 
and constructing Pade approximant with six terms in the numera- 
tor and denominator, we find x = 0.75487654 (0.00015% error). The 
method is called "the method of small delta" (see Sect. 12. 3p . 
4. We now assume the exponent as a large parameter. Consider the 
equation 

x n = 1. 



Assuming n — >• oo (a method of large J, see Sect, 
the desired solution in the form 



x = 



v 

n 



■ x\ + x 2 + . . . ) 



1/n 



(8) 

we represent 
(9) 



where 1 > x\ > x 2 > 

Substituting the ansatz (|9]) in Eq. (|8]) and taking into account that 



1/n 



1 H — In n 

n 



1 



1 + - ln(l + X! + x 2 + ...) + .. . 



one obtains in order of increasing accuracy of the formula 

m n x 

n 



In n — In In n 



n 



1/n 



(10) 



(ii) 



For n = 2 formula (fit))) gives x = 0.58871; the error compared to the 
exact solution (0.5(v / 5 - 1) w 0.618034) is 4.7%. When n = 5, we 
obtain x = 0.79715 from Eq. ([TO]) (from the numerical solution one 
obtains x = 75488; error 4.4%). For n = 5, Eq. (JlTD gives x = 0.74318 
(error 1.5%). Thus, even the first terms of the large S asymptotics give 
excellent results. 

Hence, in this case the method of large delta provides already a 
good accuracy at low orders of perturbation method. Approximation 
([TOlh ([TT]) give an example of nonpower asymptotics. 



5 



2.2 Homotopy Perturbation Method 

In recent ye ars the so-called homotopy perturbation method was popu- 
lar ( LiaoL I2QQ4} ) (term "method of artificial small parameters" is also used). 



Its essence is as follows. In the equations or boundary conditions a param- 
eter e is introduced so that for e — one obtains a BVP which admits 
a simple solution, and for e = 1 one obtains the governing BVP. Then 
the perturbation method to e is apllied and in the result ing solution was 
adopted. O f course, it i s no t new a nd was al r eady used by iGiacaglial ( 1972 



chapt. 1.6) iLiapunovl (|l893l ) and IPoincard (|l993[ ). The novelty (and, in 
my opinion, successful) is the title, emphasizing the continuous transition 
from the initial value e = to the value of e = 1 (homotopy deforma- 
tion). Let us analyse an example of homo t opy p erturbation parameter 
method using lAndrianov and Danishevs'kvvl ((20021 ) . A special feature of 



nonlinear systems with distributed parameters is the possibility of internal 
resonance between modes. That is why in many cases the neglection of 
higher modes can lead to significant errors. The following describes the 
asymptotic method of solving problems of nonlinear vibrations of systems 
with distributed parameters, allowing approximately to take into account 
all modes. The oscillations of a square membrane on a nonlinear elastic 
foundation can be written as: 

d 2 w d 2 w d 2 w a , x 

d^ + w~^- cw - £W °' (12) 

where e is the dimensionless small parameter (e <C 1). 
The BCs are as follows: 



w \x=0,L = w \y=o,L = 0. (13) 
The desired periodic solution must satisfy the periodicity conditions 

w(t) =w(t + T), (14) 

where T = ^ is the period, and ft is the natural frequency of oscillation. We 
seek the natural frequencies corresponding to these forms of fundamental 
vibration frequencies at which the linear case (e = 0) is realized by one 
half- wave in each direction x and y. We introduce the transformation of 
time: 

r = ojt. (15) 
The solution is sought in the form of expansions 

w = wo + ew\ + e 2 W2 + . . . , (16) 



6 



uj = uo + eui + £ 2 cj 2 + (17) 

Substituting ansatzes (fTBj) , (fT7|) into Eqs. ([12]) - ([T^]) and equating terms of 
equal powers , we obtain a recurrent sequence of linear BVPs: 

^ + ^--o^--o = 0, (18) 

d 2 w 1 d 2 w 1 2 d 2 w 1 d 2 w 3 , . 

^T + ^-^i^-^i = a^x-^+ti;,,, (19) 

The BCs (fT3|) and periodicity conditions (|T^|) take the form for z = 
1,2, 

^i(r) = iy*(r + 2tt). (21) 
The solution of Eq. ([T8]) is as follows: 

CO OO / \ 

^o,o = ^2 ^2 Am > n sin ( ^^-r j sin sin (~£~2/) ' ( 22 ) 

m=l m=l ^ ^ ' 



™4=o,l = °> (2°) 



where cj m)U = y7r 2 ^-^-^+c, m, n = 1, 2, 3,... and is the 
amplitude of the fundamental tone of vibrations; A m?n , m, n = 1, 2, 3..., 
(m,n) ^ (1,1) is the amplitude of the subsequent modes; u myn are the 
natural frequencies of the linear system, ujo = c^i. 

The next approximation is a result of solving the BVP ([T9]) - ([2T]) . To pre- 
vent the appearance of secular terms in the r.h.s. of Eq. ([T9]) the coefficients 
of the terms of the form 

• ( Um,n \ . firm \ . (Tin \ 

sin I —r J sin y — j^ x ) sm \~J^y) ' wi, n = 1, 2, 3, . . . 

should be equated to zero. 

These conditions lead to an infinite system of nonlinear algebraic equa- 
tions: 

~ . 000000000000 

i=l j—1 k=l 1=1 p=l s=l 

where m, n = 1, 2, 3, 

Coefficients are found by substituting ansatz ([22]) into the r.h.s. of Eqs. 
([T9]) and the relevant simplifications. System ([23]) can be solved by reduc- 
tion. However, a sufficiently large number of retention equations beging to 
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show significant computational difficulties. In addition, this approach does 
not take into account the influence of higher oscillation modes. Therefore, 
further we use the homotopy perturbation method. 

On the right side of each (m, n)-th equation of system ([23]) we introduce 
the parameter \i to those members of AijAkjA PjS , for which the following 
condition is valid: (i > m) U (k > m) U (p > m) U (j > n) U (Z > n) U (s > n). 
Thus, for ji = system ([23]) takes the "triangular" appearance, and for 
ji = 1 it returns to its original form. Next, we seek a solution in the form 
of expansions: 

wi = u/°)+^W+^V 2 ) + ..., (24) 

An,n = ^-m]n + M^-m^n + ^ ^m]n + • ■ ■ j (25) 

where m, n = 1, 2, 3, ... , (m, n) ^ (1, 1). 

We put \i = 1 . This approach allows us to keep any number of equations 
in system ([23]) . Below we limit ourselves the first two terms in the expansions 
(O, ([25]) . We analyze the solutions and note that in this problem the 
parameter c plays the role of the bifurcation parameter. In general, for 
c ^ 0, c ~ 1 , the system ([23]) admits the following solution: 



Aij, i,j = 1,2,3,..., ^ (m,n), 

27 luo 
Ul = Too ' rn,n = 1,2,3, .... 

Amplitude-frequency response is given by: 

A 2 

On,n = ^ m ,n + 0.2109375^^6: + . . . . 



Of particular interest is the case when the linear component of the restor- 
ing force is zero (c = 0) and the phenomenon of internal resonance between 
the modes of oscillations occurs. Solving system ([23]) by the described 
method, we find 

A m , n = 0, m, n = 1,2,3,..., (m, n) ^ (1, 1), (m, n) ^ (2i - 1, 2i - 1), 
i = 1, 2, 3 . . . , A 3 , 3 = -4.5662 • 10- 3 A M , A 5 , 5 = 2.1139 • 10- 5 A M , . . . , 
uji = 0.211048^^/^0. 

If the oscillations are excited by the mode (1, 1) all odd modes (3, 3), (5, 5) 
etc. are also realized. However, if the oscillations are excited by one of the 
higher modes, the result of redistribution of energy modes appear at lower 
orders until the fundamental tone (1,1). 
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2.3 Method of Small Delta 



Bender et al.l (| 19981 ) proposed an effective method of small 5, which we 



explain by examples. Let us construct a periodic solution of the following 
Cauchy problem 

x tt + x 3 = 0, (26) 
a?(0) = l, a? t (0) = 0. (27) 
We introduce a homotopy parameter 6 in Eq. ([26|h thus: 

x tt + x 1+25 = 0. (28) 

At the final expression one should put 5 = 1, but in the process of solving 
one must suppose 5 <C 1. Then 

x 26 = 1 + S\nx 2 + 0.55 2 (lnx 2 ) 2 + . . . . (29) 
We suppose the solution of Eq. ([26]) in the form 



X> fc a*, (30) 



X = 

and carry out the change of independent variable 

t=~, (31) 

where co> 2 = 1 + a\8 + a2^ 2 + . . . . 

The constants ai (i = 1,2,...) are determined in the problem solving 
process. After substituting ansatzes (|29|) - (|3T|) into Eq. (|28|) and splitting 
it with respect to 5, we obtain the following recurrent sequence of Cauchy 
problems: 

xott + x = 0, (32) 
z (0) = l, x 0r (0) = 0; (33) 

X lTT +Xi = -Xq ln(^o) - OilXOrr, (34) 

xi(0)=xi r =0; (35) 
X2tt + %2 = —xi ln(x 2 ) - 2xi - xq (\n(xl)) 2 - a 2 Xorr ~ OL\X\ TT \ (36) 

a?i(0)=xi T = 0; (37) 
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A Cauchy problem of the zero approximation ([32]) , ([33]) has the solution . 

Xq = COST. 

In a first approximation, one obtains 

x\ TT + x\ = — cost ln(cos 2 r) + ol\ cost = Lo- 

The condition of the absence of the secular terms in the solution of this 
equation can be written as follows: 

tt/2 

Lq cost dt = 0, 



which is determined by the constant a\ = 1 — 2 In 2. The expression for the 
oscillation period can be written as 

T = 27r[l + £(ln2-0.5)]. 

For S = 1, we have T = 6.8070, while the exact value is T = 7.4164 
(the error of approximate solutions is 8.2%). The solution of the Cauchy 
problems of the next approximation ([36]) , ([3?]) gives the value of the period, 
it practically coincides with the exact one (T = 7.4111). 

We now consider the wave equation 

u u = u xx , (38) 

with nonlinear BCs: 

u(0,t) = 0, (39) 

u x (l,t)+u(l,t)+u 3 (l,t) = 0. (40) 

We introduce the parameter S in Eq. ([^0]) as follows: 

1^(1, t) + tx(l, t) + u 1+25 (l, t) = 0. (41) 

In the final expression should we put 5 = 1, but in the process of decision 
we assume S <C 1. Then 



u 3 = u 1+2S 



1 + S In u 2 + — (In u 2 



(42) 



Suppose further the solution of Eq. ([38]) in the form 

oo 

u = J2^u k . (43) 



fc=0 
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After substituting ansatzes (|43]), (@TJ) into Eqs. (|38]h (l39l) . (|4B and the 
splitting of the parameter 5 we obtain the following recurrent sequence of 
BVPs: 



at x = 0, wo = 0; 

at x = 1, iiox + 2iio = 0; 
l 



p=0 

at x = 0, ?ii = 0; 
at x = 1, uia; + 2^i = —^o In 

2 

^Orr — V>0xx ^ ^ ^i—p^prri 

at x = 0, = 0; 
at x = 1, + 2^2 = — ui ln^Q — 2^i — 0.5^o (ln^o) 



(44) 
(45) 
(46) 

(47) 

(48) 
(49) 

(50) 

(51) 
(52) 



where ao = 0. 

The solution of the BVP of the zero approximation can be 

written as 

uq = Asm(ujQx) sin(co>oT), (53) 



where the frequency ujq is determined from the transcendental equation 



coo = 2 tanco>o- 
The first few nonzero values of uo are given in Table [TJ 



(54) 













2.289 


5.087 


8.096 


11.173 


14.276 












, ,(10) 


17.393 


20.518 


23.646 


26.778 


29.912 



Table 1. The first few roots of the transcendental Eq. (j5l]). 



11 



When k — » oo we have the asymptotics: uj^ — » 0.57r(2/c + 1). 
The BVP problem of the first approximation are as follows: 

uixx — ui TT — ol\Au\ sin(cjo^) sui(co>ot); (55) 
at x = 0, u\ = 0; (56) 
at x = 1, uiaj + 2^i = Ai sin(cJoT) [ln(A 2 sin 2 cJo) + lnsin 2 (co>oT)] ; (57) 

where A\ — — AsincJo- 

The particular solution of Eq. (|55|) satisfying the BC ([56]) has the form 

= — 0.5aiAo;o^cos(a;o^) sin(cJoT). (58) 

We choose the constant ol\ in such a way that it compensates the secular 
term on the r.h.s. of Eq. (j57|) 

2i?i 

ai = 



cjo(6 + cJq) 



where i?i = ln(0.25eA 2 sin 2 cj ). 

Non-secular harmonics on the r.h.s. of Eq. (|57j) give the solution 

u^ 2) = 4Ai T k sin(o;Qfex) sin(cj fer) — - , (59) 

k=2 

where Tk = l/[kuo cos(kuo) + 2sin(fcu;o)]- 

The complete solution of the first approximation has the form 

(1) , (2) 

Assuming S = 1, we obtain the solution of Eqs. (I38l)-(l40l)_. 

ion (|Boe 



Let us now consider the Schrodinger equation ([Boettcher and Bender , 

1990h 

*xx - x 2N ^ + E9 = 0, (60) 

tf(±oo)=0. (61) 

Here \£ is the wave function; E 1 is the energy,it plays the role of eigen- 
value. It is shown that the eigenvalue problem (1501) 2 (15T1) has a di screte 
countable spectrum E n , n = 0,1,2,... ([Slepvan and Yakovlevl |l980). For 
N = 2 the eigenvalue problem ([60]) , (|5T]) has an exact solution. Now let TV 
differ little from 2, 

^ - x 2+25 ^ + ^ = 0. (62) 
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Now we use the expansion 

x 26 = l + 5\n(x 2 ) + ..., 

and we will search for the eigenfunction \I/ and the eigenvalue E in the form 
of the following series: 

^ = ^ + ^ + £ 2 ^ 2 + ..., (63) 

E = E + SE + S 2 E 2 + . . . . (64) 

As a result, after the asymptotic splitting we obtain a recursive sequence of 
eigenvalue problems 

*Oxx-x 2 *o + £o^o = 0, (65) 
^i xx - £ 2 ^i + £0^1 + £1^0 = a 2 * ln(x 2 ), (66) 



l^il-^O at \x\->oo, i = 1,2,3,.... (67) 
The solution of the eigenvalue problem ([65]) , ([67]) has the form 

4 n) = 2n + 1, = e-* 2 / 2 tf n (;r), n = 1, 2, 3, . . . , 



where H n (x) is the Struve function ([Abramowitz and StegunlTl965i Sect. 12). 



From the eigenvalue problem ([66]) . (|67|) we find 



gin) _ -°° 



jf x 2 e x<2 H 2 (x) \n(x 2 ) dx 



v/tt 2 n n! 

For n = one obtains Hq(x) = 1 and 



00 

/ 



x 2 In k"" 2 dx = ^ (2 - 2 In 2 - C) , 
8 



where C = 0.577215 ... is the Euler constant. Hence 

4 0) = 1 + 4(2 -21n2- + . (68) 
16 
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2.4 Method of Large Delta 

An alternative method of small delta is the method of large delta, which 
we demonstrate in an example of a nonlinear equation 

x tt +x n = 0, n = 3,5,7,... . (69) 

This equat ion can be integr ated with the functions Cs and Sn, introdu ced by 
Liapu nov iLiapunovl (|l893[ ) (inversions of incomplete beta functions, ISenikl 
(1969)). Note that much later the same (up to normalization) fun ction have 
been p roposed by Rosenberg, who called them Ateb-functions iRosenbergl 
(1963). However, working with these objects is inconvenient, and therefore 
the problem of the approximate analytical solution of Eq. ([69]) in elementary 
functions arises. We construct asymptotics of periodic solutions of Eq. ([69]) 
at n — » oo. Let the initial conditions for Eq. ([69]) be 

x(0) = 0, x(0) = 1. (70) 

The first integral of the Cauchy problem ([69]) , ([TO]) can be written as: 

V dt J n + 1 v ' 

The replacement of the x = A -A / 2 , A = 2/(n + 1) , and integration gives us 
a solution in the implicit form 

o<£<i 



/ 



V^W^' 

After replacing ^ = sin A this implicit solution is transformed into an 
expression that contains a small parameter in the exponent of the integrand: 

o<e<7v/2 

X x/2 t = X I sin- 1+A 6>d<9. 



o 



We now consider the integrand separately: 



sm- 1 + A ^ = ^- 1+A f — > -«" 1+A 
sin^ 



Aln- 



sm V sin t 



Expanding this function into a Maclaurin series, one obtains 

sin- 1+A e = e~ l+x + + * * ' + °( A )- 
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The first term of this expression makes the main contribution of this 
expression, so in the first approximation we can write 



i.e. 



A l/2 t l/A - 



In the original variables one obtains 

x « A" A /2 sin A (a 1 / 2 ^) (72) 
The solution (f72l should be used on a quarter-period Twith 



(2AV2) 



(73) 



Let us analyze the solution ([72]) , ([73]) . At n = 1 one obtains the exact 
values x = sint, T = 27r , for n — ^ 00 one obtains T — » 4. Expanding the 
r.h.s. of the Eq. ([?2]) in a s eries of t and res tricting it to the first term, we 
obtain nonsmooth solution ( PilipchukL 2010[ ). We estimate the error of the 
solution ([73]) . For this we use the expression 



t/2 



-1+A 



dO = 0.5A£(0.5A,0.5) = A u 



(74) 



where B(. ..,...) is the beta function ( Abramowitz and Stegunill965 , Sect. 6). 
The approximate value of the integral on the l.h.s. of Eq. ([72]) is calculated 
as follows: A2 = (tt/2) a . Numerical comparison of the values Ai, A2 and 
error estimate A is given in Table [2j 



n 


1 


3 


5 




00 


A! 




1.30 


1.20 




1 


A 2 




1.25 


1.16 




1 


A,% 





5 


3 




- 



Table 2. Comparison of exact and approximate solutions. 



Thus, the first approximation of the asymptotics for n —> 00 already 
gives quite acceptable accuracy for practical purposes, even for not very 
large values of n. Note that the expression (l72l) give s an a pproximation 
of incomplete beta function ([Abramowitz and Stegunl 
n = 1 (sinus function) to n = 00 (linear function). 



19651 Sect. 6) from 
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2.5 Application of Distributions 

Asymptotic methods are based, generally speaking, on the use of Taylor 
series. In this connection the question arises: what to do with the functions 
of the form exp(— which c annot be exp anded in a Taylor series in 
e — » 0, using smooth functions (|PeierlsjJl979h. The way out lie s in the 



transition to the following distribution ( Estrada and Kanwali I2OO2I ) 



H(x)exp(-e- 1 x) = ^(-l) n £ n+1 J^(x), (75) 



n=0 



where S(x) is the Dirac delta function, representing the derivative of the 
Heavisyde function H(x) in the theory of distribution; 5( n \x), n = 1, 2, . . . 
are the derivatives of the delta function. 

We show how formal the formula ([75]) can be obtained. Applying the 
Laplace transform to function exp(— e~ 1 x), one obtains: 



00 

/ 



exp(— e x x) dx — 



ep+1 



Expanding the r.h.s. of this equation in a Maclaurin series of e, and then 
calculating inverse transform term by term, one obtains expansion ([75]) . 
Thus, we again use the Taylor series, but now in the dual space. 

Here is another interesting feature of the approach using distributions: 
a sing ular perturbated pr oblem can be regarded as a regular perturbated 
one ([Estrada and KanwalL 20021 ). Suppose, for example, 



ey' + y = x > 0, y = 1 at x = 0. 

This is singular perturbated problem: for £ = one obtains a smooth 
solution y = 0, which does not satisfy the given initial condition. However, 
one can seek a solution in the form of a nonsmooth function. Namely, 
assuming z(x) = H(x)y(x), one obtains 

ez' = -z + e5(x). (76) 

The solution of Eq. ([76]) is sought in the form of expansions 

00 

n=0 

As a result, one obtains 

z = 0, Zl = 5(x), z n+1 = (-l) n 5^ n \ n = l,2,.... (77) 



16 



Note that the expressions ([77]) can go to the smooth functions. To do this, 
it is possible to apply the Laplace transform, then the Pade approximants 
in the dual spase and then to calculate inverse Laplace transform. 

We show other application of the asymptotic method using distribution. 
Consider the equation of the membrane, reinforced with fibers of the small, 
but finite width e. The governing PDE is: 

[l + 2e&o(y)]u xx +Uyy = 0, (78) 



where $ (y) = ]T [H(y + kb - e) + H(y - kb + e)]. 

k=—oo 

Let us expand th e function <I>o (y) in a series of e. Applying the two-sided 



Laplace transform ( van der Pol and Bremmerl Il987 ), one obtains: 



oo 

*(p,e) = J e- plyl $(y,e) dy 



Expanding the function 3>(p, s) in a series of e and performing the inverse 
Laplace transform, we obtain 

*o(y) = 2 £ $(y) + 2 £ Y, e " $(n) (y), (79) 

fc=l,3,5,... 



where $>(y) = XI ^(v ~ kb). 

k= — oo 

Now let us consider the solution of the Eq. ([78]) in the form 

u = + sui + s 2 u 2 + (80) 

Substituting ansatzes ([79]) . ([80]) into Eq. ([75]) and splitting the resulting 
equation with respect to e, we arrive at the recursive sequence of BVPs 

[l + 2e$(2/)]u 0a;a; + uo 3/2/ = 0, (81) 
[l + 2£$(y)]u lxx + u lyy = -eu 0xx $ y (y), (82) 



Thus, in the zero approximation, we obtain the problem with one-dimensional 
fibers ([8T]h the influence of the width of the filaments is taken into account 
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in the following approximations ([82)). 



3 Summation of Asymptotic Series 

3.1 Analysis of Power Series 

We suppose that one obtains the following series as a result of an asymp- 
totic study: 

/(e)~5^C n e n for £ ^0. (83) 

n=0 

As it is known, the radius of convergence eo of series ([83]) is determined by 
the distance to the nearest singularity of the function f(e) on the complex 
plane and can be found using the Cauchy-Hadamard formula: 

-=TnrT|C„| 1/n . 

If the nearest singularity lies on the positive real axis, then the coefficients 
C n have the usually one and the same algebraic sign, for example, 

1 ~l+£ + £ 2 +£ 3 + .... 



1 - 

If the nearest singularity is located on the negative axis, the algebraic signs 
of the coefficients C n are usually alternated, for example, 

1 ~ 1 - e + e 2 - e 3 + . . . . 



1 + 

The pattern of signs is usually set pretty quickly. If there are several features 
of the same radius, that could happen to a real function with complex 
singularities necessarily occurring in complex conjugate pairs, then the rule 
of alternation of signs may be more complex, such as 



Here we have a pattern of signs -H - . To define sp may it is useful t o apply 

the so-called Domb-Sykes plot (iHinchi Il99l[ Ivan DvkeL Il97l Il975al lbh. Let 
the function / have one of the nearest singularities at a point e = ±£o with 
an index of a, i.e. 



(e ±e) a for a ^0,1, 2,..., 

(s ± e) a ln(e ± e) for a = 0, 1, 2, . . . , 
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then we get 

C n ^ l_ L _ l + a 

C n -i £o \ n 

Constructing a graph of C n /C n -\ an the vertical axis and 1/n on the 
horizontal axis, one obtains the radius of convergence (as the reciprocal 
of the intercepts on the axis C n /C n _i), and then, knowing the slope, the 
required singularity. Figure 1 shows the numerical results for the function 



f(e) = e(l + e)(l + 2e)- 



-1/2 



9 3 „ 3 4 27 ~ 51 fi 191 7 359 R 

^ p _ p2 j £ 3 £ 4 j 5 p6 _, £ 7 £ 8 , 

+ 2 2 + 8 8 + 16 16 £ 

starting with n = 7 points arranged in a linear relationship. 



















> o 
























O j 








1 









Figure 1. The Domb-Sykes plot for f(e) = e(l + e)(l + 2s)- 1 / 2 . 



If So or <^ are known from physical considerations, they can be used 
for the construction of the Domb-Sykes plot. If several singularities have 
the same convergence radius, so that the signs of the coefficients oscillate, 

1/2 

you can try to construct a dependence on the value (C n /C n _i) . If the 
radius of convergence tends to infinity and C n /C n -\ ~ k/n, then the ana- 
lyzed function has a factor exp(fcs), when C n /C n -i ~ k/n 1 ^ it has a factor 
exp(£ p ). If the radius of convergence tends to zero, then the analyzed func- 
tion has an essential singularity and asymptotic expansion diverge. If the 
coefficients behave like C n -\/C n ~ l/(kn) then we can write C n ~ Ck n n\, 
where C is a constant. 



19 



Knowledge of the singular solutions can eliminate them from the pertur- 
bation series and thus significantly improve its convergence. We describe 
some techniques for removing singularities. If the singularity lies on the 
positive real axis, then it often means that the function f(e) is multivalued 
and that there is a maximum attainable point e = Sq. Then the inverse of 
the original function e = e(f) can be single valued. For example, consider 
the function 

f(e) = arcsine = e + ^ 3 + + -^e 7 + . . . , (85) 

the inverse of this function is 

£ ~ f ~ -f + -20 f — f + . . . . (86) 

Numerical results are shown in Figure 2, where the solid line denotes the 
function arcsine, the dotted and dashed line shows the n-term expansions 
(|85j) and k-terms expansions ([86]) for different numbers of terms. It is evident 
that the expansion ([86]) allows a good description of the second branch of 
the original function. 




Figure 2. Application of the inversion of a power series. 



If 

f~A(e -s) a for e e , < a < 1, 
the transition to the function J 1 / 2 removes the singularity. 
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Consider the following example: 

e/2 r FT" , 1 7 2 41 3 367 4 4849 5 

fie) = e-^VT+Ye -l + - 2 e--/ + T / - _ e < + _*» + ... . (87) 

The radius of convergence of this expansion is equal to 1/2, while the 
radius of convergence of functions 

is infinite. 

Numerical results are shown in Figure 3, where the solid line denotes 
the function f(e) = e~ e l 2 \J\ + 2e, the dotted and dashed lines shows the 
n-term expansions (|8TJ) and the square roots of /c-term expansions ([88]) . 

u 




0-5" I , I . — —J 1 — . i-^ — - — . — ■ — L-#- 

0,0 as to is £ 

Figure 3. An example of taking the root. 

In addition, knowing the singularity, one can construct a new function 
/m(^) (multiplicative extraction rule) 

/(e) = ( £o - e) a f M (e), 

or /a(^) (additive extraction rule) 

f(e) = A(e - eTf A (e). 

The functions /m(^) an d /a(^) should not contain singularities at £q. 
In many cases, one can effectively use the conformal transformation of the 
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series, a fairly complete catalog of which is given in (|Kublanovskava.l Il953h . 
In particu l ar, it sometimes t urns ou t to be a successful Euler transformation 
(iBellmanl . Il964i : Ivan Dyke] , Il975a lbl), based on the introduction of a new 
variable 

Recast the function / in terms of e , / ~ ^2d n s n , has the singularity 
pushed out at the point e — oo. For example, the function (|8^|) is singular 
at the e = —1/2, which can be eliminated with the Euler transformation 
£ = £/(! + 2e). The expansion of the function (|51]) in terms of e is 



1 



/(g) „ i + g + - — e d - — 



31 o 895 , 22591 , 



(90) 



2 8 48 384 3840 
Some numerical results are shown in Figure 4, where the dotted and dashed 
line shows the n-term expansions (|8^|) and /c-terms in the expansion ([90]) . 



/ 



7.0 



0.5 





\ **' . ' / 

— 






+ « 
* 

i, i — i — i — : 





K*6 



Figure 4. Illustration of Euler transformation. 
A natural generalization of Euler transformation is 

£ 



e = 



(l-e/e ) a ' 



where a is the real number. 



3.2 Pade Approximants and Continued Fractions 

"The coefficients of the Taylor series in the aggregate have a lot more 
information about the values of features than its partial sums. It is only 
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necessary to be able to retriev e it, and some of the way s to do this is to con- 



struct a Fade approximant" ([Vinogradov et all Il987l ). Pade approximants 



(PA) allow us to implement among the most salient natural transformation 
of power series in a fractional ra tional function. Let us define a PA following 
iBaker and Graves-Morris! ( 19961 ). Suppose we have the power series 



oo 



f(e) = ^c z e\ (91) 

= 1 



Its PA can be written as the expression 

a + ai£H Va n e n 

f ^ {£) = l + 6ie + ... + 6m£ m> ( 92 ) 

whose coefficients are determined from the condition 

• ■+b m £ rn )(co+c 1 £+c 2 s 2 + . . . ) = a +ai£+- • ■+a n s n +0(e m + n + 1 ). 

(93) 

Equating the coefficients of the same powers e, one obtains a system of 
linear algebraic equations 



+ c n+ i 0; 

b m £n-rn+2 + ^m-lC n -m+3 + C-n+2 = 0; 
bm^n H~ bffi— iC n _|_i -|- C n _|_ m — 0, 



(94) 



where Cj = for j < 0. 

The coefficients a$ can now be obtained from the Eqs. (|93j) by comparing 
the coefficients of the powers e: 



a = c ; 

ai = ci + 6ic ; 



i=l 



(95) 



where p = min(n, m). Eqs. ([93|) . ([M|) are called Pade equations. In the case 
where the system ([94]) is solvable, one can obtain the Pade coefficients of 
the numerator and denominator of the PA. Functions f[ n / m ] (e) at different 
values of n and m form a set, which is usually written in the form of a table, 
called the Pade table (Table 3): 

The terms of the first row of the Pade table correspond to the finite sums 
of the Maclaurin series. In case of n = m one obtains the diagonal PA, the 
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n 

m ^\ 





1 


2 







ijo/oiO) 


/ri/oi (^) 


/[2/01 O) 




1 


/fo/il(^) 


/fi/nte) 


/[2/n (e) 




2 


/f0/2l(e) 


/[1/21OO 


/[2/21 0) 















Table 3. Pade table. 



most common practice. Note that the Pade table can have gaps for those 
indices n, m, for which the PA does not exist. We note some propertie s of 
the PA (|Apresvani 1 19791 : iBaker and Graves-Morrisl Il996l : ISuetinl . [20021 ): 



1. If the PA at the chosen m and n exists, then it is unique. 

2. If the PA sequence converges to some function, the roots of its denom- 
inator tend to the poles of the function. This allows for a sufficiently 
large number of terms to determine the pole, and then perform an 
analytical continuation. 

3. The PA has meromorphic continuation of a given power series func- 
tions. 

4. The PA on the inverse function is treated the PA function inverse 
itself. This property is called duality and more exactly formulated as 
follows. Let 

q(e) = f- 1 (e) and /(0) + 0, then q [n/m] (e) = f^ /m] (e) , (96) 

provided that one of these approximations there. 

5. Diagonal PA are invariant under fractional linear transformations of 
the argument. Suppose that the function is given by their expansion 
(|9T|) . Consider the linear fractional transformation that preserves the 
origin W = -f^, an d the function q(W) = f(s). Then q[ n / n ], pro- 
vided that one of these approximations exist. In particular, the diag- 
onal PA is invariant concerning Euler transformation ([89]) . 

6. Diagonal PA are invariant under fractional linear transformations of 
functions. Let us analyse a function (|9T|) . Let 



If c + df{e) ^0, then 

Q[n/n] (e) 



df(e)' 

• bf [n / n ] (e) 
■ df[ n / n ] (e) ' 
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provided that there is f[ n / n ] (e). Because of this property infinite values 
of PA can be considered on a par with the end. 

7. The PA can get the upper and lower bounds for f[ n / n ](s). For the 
diagonal PA one has the estimate 

/[n/n-l](e) < f[n/n]{e) < /[n/n+l](^)« (97) 

Typically, this estimate is valid for the function itself, i.e. f[ n / n ](s) 
in Eq. ([97|) can be replaced by f(e) . 

8. Diagonal and close to them a sequence of P A often possess the prop- 
erty of autocorrection ( Lit vinovl 1 1 994 [2003 ) . It consists of the follow- 



ing. To determine the coefficients of the numerator and denominator 
of PA have to solve systems of linear algebraic equations. This is 
an ill-posed procedure, so the coefficients of PA can be determined 
with large errors. However, these errors are in a certain sense of self- 
consistency, the PA can approximate the searching function with a 
higher accucary. This is a radical difference between the PA and the 
Taylor series. 

Autocorrection property is verified for a number of special functions. At 
the same time, even for elliptic functions the so-called Froissart doublets 
phenomenon arises, consisting of closely spaced zeros and poles to each 
other (but different and obviously irreducible) in the PA. This phenomenon 
i s not of num erical nature, but due to the nature of the elliptic function 
( SuetinLl2002l ). Thus, in general, having no information about the location 



of the poles of the PA, but relying solely on the PA (computed exactly as 
you wish), we can not say that you have found a good approximation for 
the approximated function. 

To overcome these d efects several methods are suggeste d, in particular, 



the smoothing method (|Beckermann and Kaliaguind 119971 ). Its essence is 



that instead of the usual-term diagonal PA for complex functions f[ n / n ] (e) = 
Pn(z) / Qn(z) the following expression is used 

, ( \ — ^n(g)Pn(g) + gn-l(g)Pn-l(g) 
J[n/n]\ e ) -r^ , . -r^ , . 

Here / denotes complex conjugation of /. 

Now consider the question: in what sense can the available mathemati- 
cal results on the convergence of the PA facilitate the solution of practical 
problems? Gonchar's theorem (Gonchar, 1986) states: if none of the di- 
agonal PA f[ n / n ](e) has poles in the circle of radius R, then the sequence 
f[ n / n ](s) is uniformly convergent in the circle to the original function f(s). 
Moreover, the absence of poles of the sequence of the f[ n / n ](s) in a circle 
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of radius R must be original and confirm convergence of the Taylor series 
in the circle. Since for the diagonal PA invariant under fractional linear 
maps we have e — >• ^q^, the theorem is true for any open circle contain- 
ing the point of decomposition, and for any area, which is the union of 
these circles. A significant drawback in practice is the need to check a ll 
diagonal PA. For example, the following theorem holds ( VvatchinL Il982[ ): 



Suppose the sequence of diagonal Pade approximants of the function w(e), 
which is holomorphic in the unit disc and has no poles outside this circle. 
Then this sequence converges uniformly to w(e) in the disc \z\ < 7*0, where 
0.583i? < r < 0.584R. 

How can we use these results? Suppose that there are a few terms of the 
perturbation series and someone wants to estimate its radius of convergence 
R. Consider the interval [0, £o]> where the truncated perturbation series and 
the diagonal PA of the maximal possible order differ by no more than 5% 
(adopted in the engineering accuracy of the calculations). If none of the 
previous diagonal PA does not have poles in a circle of radius Sq : then it is 
a high level of confidence to assert that R> e$. 

The procedure of constructing the PA is much less labor intensive than 
the construction of higher approximations of the perturbation theory. The 
PA is not limited to power series, but to the series of orthogonal polynomi- 
als. PA is locally the best rational approximation of a given power series. 
They are constructed directly on its odds and allow the efficient analytic 
continuation of the series outside its circle of convergence, and their poles 
in a certain sense localize the singular points (including the poles and their 
multiplicities) of the continuation function at the corresponding region of 
convergence and on its boundary. This PA fundamentally different from 
rational approximations to (fully or partially) fixed poles, including those 
from the polynomial approximation, in which case all the poles are fixed 
in one, infinity, the point. Currently, the PA method is one of the most 
promising non-linear methods of summation of power series and the local- 
ization of its singular points. Including the reason why the theory of the 
PA turned into a completely independent section of approximation theory, 
and these approximations have found a variety of applications both directly 
in the theory of rational approximations, and in perturbation theory. Thus, 
the main advantages of PA compared with the Taylor series as follows: 

1. Typically, the rate of convergence of rational approximations greatly 
exceeds the rate of convergence of polynomial approximation. For ex- 
ample, the function e £ in the circle of convergence approximated by 
rational polynomials P n (s)/Q n (e) in 4 n times better than an algebraic 
polynomial of degree 2n. More tangible it is property for functions 
of limited smoothness. Thus, the function \e\ on the interval [—1,1] 
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can not be approximated by algebraic polynomials, so that the or- 
der of approximation was better than 1/n, where n is the degree of 
polynomial. PA gives the rate of convergence ~ exp ( K —^/2n). 

2. Typically, the radius of convergence of rational approximation is a 
large compared with power series. Thus, for the function arctan(x) 
Taylor polynomials converge only if \e\ < 1, and AP - everywhere in 
C \ ((—zoo, —i] U [i, ioo)). 

3. PA can establish the position of singularities of the function. 
Similarly, the PA method is a method of continued fractions (Uones and Throni 



1980). There are several types of continued fractions. The regular C 
fraction has the form of an infinite sequence, in which 7V-th term can be 
written as follows 

f N (e) = a + \- e . (98) 



c 2 s 



1+ 

cn+i£ 



I + C N £ 

The coefficients q are obtained after the decomposition of expression 
([98|) in a Maclaurin series and then equating the coefficients of equal powers 
of e. When a = one obtains the fraction of Stieltjes or S- fraction. For the 
function of Stieltjes 

oo 



the coefficients of expansion ([98]) have the form: a = 0, Co = 1, C2 n -i = 
C2n = n. n > 1. 

Description of the so-called J-, T-, P-, R-, g- fractions, algorithms for 
their construction and the range of applicability are described in detail in 
Jones and Throni (|l980l ). 



Conti nued fractions are a specia l case of continuous functional approx- 
imation (Bender and Orszad . Il978h . This is the sequence in which the 



(n + l)-th term c n (e) has the form n-th iteration of a function F(e). For 
the Taylor series one obtains F(e) = 1 + e, for the continuous fraction 
F( £ ) — i+p ^ F( £ ) = exp(e) one obtains a continuous exponential approx- 
imation 

c n (e) = a exp{aieexp [a 2 e . . .exp (a n e)]} , 
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for F(e) = vTTi 



c n (e) = a Y 1 + aiey 1 + a 2 £\/ 1 + . . . a n _i£Vl + a n £, 
for = In 1 + e 

c n (e) = a ln{ai£:ln[a2£:...ln(a n £:)]}. 

In some cases, such approximations can converge significantly faster than 
power series. 

As an example, we note the solution of the transcendental equation 

x = e In x 



for large values ([ArgatovL 12004 §3.4.9)) A: 



xo=ehie] x\ =e Inline); X2 = s\n[s\n(s\ns)}; 

4 Some Applications of Pade Approximants 

4.1 Accelerating Convergence of Iterative Processes 

The efficiency of PA or other methods of summation depends largely on 
the availability of higher approximations of the asymptotic process. Some- 
times they can be obtained by using computer algorithms (Millerl 120061 ), 



but in general it remains an open question. Iterative methods are signifi- 
cantly easier to implement. As a result of an iterative procedure a sequence 
of S n is obtained. Suppose that it converges and has the limit value. We 
introduce the parameter a obtained by the ratio 

i • S'n+l — S n 

a = hm 



n— )>c>o 



S n s 



It's called superlinear convergence, if a = 0, a linear for a < 1 and 
logarithmic at a = 1. The biggest issues are, of course, logarithmically 
convergent sequences - alas, widespread practice. Very often linearly con- 
vergent sequences are also a problem. Therefore, it is often necessary to 
improve the convergence. One method of improving the convergence is to 
move to a new sequence T n with the aid of a transformation so that 

ii m T :~ s : = 0. 
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In such cases we say that the sequence T n converges faster than sequence 
S n . There are linear and nonlinear methods to improve convergence. Linear 
methods are described by formulas 



0,1,2,... 



where the coefficients a n i do not depend on the terms of the sequence S n 
constant. 

Since linear methods improve the convergence of a restricted class of 
sequences, currently nonli near methods the most popular ones. Among 
them thr Aitken method (Bak er and Graves- Morrisi 19961 ) stands out for 
its easiness, which described by the formula 



(Sn+l — S n ) (S n ~ S n -i) 



s, 



n+1 



2S n 



0,1,2,... . 



(99) 



The Aitken method accelerates the convergence of all linear and many 
of logarithmically convergent sequences. It is very easy to calculate, and 
in some cases it can be applied iteratively. A natural genera l izatio n of the 
Aitken transformation is the Shanks transformation (|ShanksL 0.955) 



rpsh 



D 



(i) 

kp 



D 



(i) : 

kp 



(100) 



where 



D 



(i) 

kp 



D 



(2) 
kp 



Sp—k 

ASp-k 
1 

AS p -k 
ASp-i 



Sp-k+i 



Sp 
ASp 

ASp+k-i 
1 

ASp 

ASp+k-i 



ASk = Sk+i — Sk 

Eq. ()100p is called the Shanks transformation of the order k of the 
sequences to the sequence Sk- For k = 1 one obtains the Aitken transform 
([99]) . Shanks method requires the calculation of determinants, which is not 
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always easy. One can use also Wynn algorithm, which described by the 
formulas: 



rp(n) _ n rp{n) q rp{n) _ rpin+l) ± /ini\ 

i -l ~ U ' i -^n, + ™(n+l) _ ^(n) ' V iUi J 

The Wynn algorithm is related to the transformation of Shanks (|100p in the 
following way: 

rp(n) _ rp(sh) ( q \ rr^n) _ ' ^ A C \ 

The Wynn algorithm is a quadratic convergent method for solving sys- 



tems of nonlinear equations ( Baker and Graves- Morris! 19961 ). There are 



many other techniques for accelerating of sequences' convergence. One can 
use them consistently, for example, to convert the original sequence into a 
linearly convergent one, and then apply the method of Aitken. One can also 
use different methods to improve conv ergence, at each stage by comparing 
the obtained results (|Brezinskil 120001 ). All the described methods have a 
close relationship with the PA. The Aitken method corresponds to the PA 
[n/1] , the Shanks method to the PA \p/k], for the method of Wynn one 
obtains = [n + k/k]. 

4.2 Removing Singularities and Reducing the Gibbs Effect 

Consider the problem of uniform plane flow of an incompressible inviscid 
fluid pasting a thin elliptic airfoil (\x\ < 1, \y\ < g, e *C ! )• The expression 
for the relative velocity q* of the flow is such ( van DvkeUl975bl §4.4): 



V ^/l-x 2 (l + s 2 ) v J 

where V is the free-stream speed. 

The splitting of the r.h.s. of Eq. (j!02[) in a series of e can be expressed 

as: 

1 J2 1 ™2 

„2 x 1 „3 



g .( x , e ) = l + e __ e i____ ed __ + ... . (103) 

This expression diverges at x — 1. No wonder it is not: the expansion 
(|103p is obtained as a result of the limiting process lim q(x,e), x > 1, and 

to get the value of g(l,g), it is necessary to perform the limit as lim q(x,e) 

x—*l 

for e > 0. Divergence of series (|103p when x = 1 indicates that the limit 
processes cannot be interchanged. Now let us apply PA to the r.h.s. of Eq. 
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(|103j) and then pass to the limit x —> 1. After trying various options, we 
conclude that the best result is given by the PA 

^)= "7 2)(1 2 +£) . (104) 
1 — or 

Numerical results for £ = 0.5 are shown in Figure 5, where the dashed 
line denotes the solution (|103[L curves 1 and 2 the exact solution (|104|) 
and the PA (jl04[) . It is seen that the use of PA significantly improves the 
accuracy of the approximate solution. 




PA ca n be also successfully appl i ed for the suppression of the Gibbs phe- 
( Becker mann et all 120081 : iBrezinskii l2000t iDriscoll and Fornberg, 



nomenon 



200ll ; [Nemeth and ParisL 1985( 1 . Consider, for example, the function sign x: 

-1, -7r < x < 0, 

1, < X < 7T. 



sign x 



Its Fourier series expansion has the form 



sign x 



3=0 



sin(2j + l)x 

2jn 



(105) 
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Direct summation of series (|105p leads to the Gibbs effect in the neigh- 
borhood of x — 0, while the defect of convergence reaches 18%, i.e. instead 
of 1 one obtains the value of 1.1789797. . .. Diagonal PA for series (|105j) 
can be written as follows: 



where 



S2t 



Sign X [AT/AT] 



N-l/2 

g2j+isin((2j + l)x) 
j=Q 

N/2 

1 + £ s< 2j cos(2jx) 



<Z2i+i = -(2j + l) 

7T 



1 



S2* 



[AT/2] 

£ (2j + l)2-(2i) 2 



= 2(-l)< 



(106) 



(2j + l)^ 

(AH) 4 (27V + 2i)!(2iV-2z)! 
(TV - 1)!(7V + 1)!(7V - 2i)\(N + 2z)![(2JV)!] 2 ' 



Numeric al studies show that the Gibbs effect for PA (|106[) does not 
exceed 2% (jNemeth and Paris! . fl985l ) . 



4.3 Localized Solutions 

We consider the stationary Schrodinger equation 

V 2 (x,y) - u(x,y) + u 3 (x,y) = 0. 



(107) 



We seek the real, localized axisymmetric solutions of the Eq. (|107j) . In polar 
coordinates (£, 0) we construct a solution which does not depend on 
As a result, we obtain the BVP 



v"(0 + |^(0-v(0 + v 3 (0 = o, 

<p(x) = 0, 
lim ip(0 = 0. 



(108) 

(109) 
(110) 



BVP (|108|) - (|110p can be regarded as an eigenvalue problem, and the role 
of an eigenvalue is an unknown quantity A = cp(0). This problem plays an 
important role in nonlinear optics, quantum field theory, the theory of mag- 
netic media. As shown in (|Molotkov and VakulenkoL Il988l p.12-16), BVP 
(|108p - (|110p has a countable set of "eigenvalues" A n , the solution ip(£,A n ) 
has exactly n zeros, and the solution ip(£,Ao) has no zeros and decreases 
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monotonically on £. That is the last solution, which is most interesting 
from the standpoint of physical applications, and we will concentrate on 
obtaining it. 

The problem of computing the decaying solutions of BVP (|108)) - (|110p 
is identical to the problem of computing homoclinic orbits in the 3D phase 
space for the nonlinear oscillator, or equivalent ly, for computing the initial 
conditions for these orbits. 

Since the sought solutions are expected to be analytical functions of £, 
they can be expressed in Maclaurin series about £ = 0: 

CO 

^)=A +^c 2J e j . (in) 

Substituting ansatz ([Hip into Eq. (|108[L producing a s plitting of the 



powers of the £ and solving the relevant equations, one obtains (|Emachi et al 



19971 ): 



C 2 = 0.25A (1 - Al)- C 4 = 0.25C 2 C; C 6 = -f C 8 = 

^(CC 6 - 6A C 2 C 4 - Cf); 

C w = -0.01(CC 8 + QA C 2 Ce + 3A C% + 3C 2 2 C 4 ); 
C12 = -^(-CC W + 6A C 2 Cs + 6A C 4 C 6 + 3C 2 c 2 4 ), 

where C = 1 - 3A%. 

Then we construct PA for the truncated series 

Ao + E a 2j e j 

<p(Z) = j=5 (H2) 

i + e b 2J e k 

k=l 

All coefficients in Eq. (|112p can be parameterized in terms of Ao, <22j = 
a,2j(Ao), b 2 j = b2j(Ao) and the PA becomes a one-parameter family of 
analytical approximations of the searching solution. Then we compute the 
value of Aq for which PA (|112|) decays to zero as £ tends to infinity. It gives 
us condition 

02j(Ao) =0, 6 2j -(A ) ^0, (113) 

One can compute the PA (|112p , then imposed the condition (|113p and 
obtained the following convergent values of Aq for varying orders 27V: 
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N 


1 


2 


3 


4 


A 


±Vs 


±2.20701 


±2.21121 


±2.21200 



The numerical solution gives Aq « ±2.206208, the difference between 
numerical and analytical solutions for N = 4 is 0.26%. 

4.4 Hermit e-Pade Approximations and Bifurcation Problem 

PA can successfully work with functions having poles. However, it often 
becomes necessary to explore functions with branch points, and construct 
all their branches. In t hat c ase, one can use Hermite-Pade approximations 
( Drazin and TourignvL Il996h . Suppose it comes to a function with the ex- 
pansion 

CO 

/(e) = 5> n e n , (114) 

n=l 

and we managed to find the first few coefficients of this series 

N 

fN(e) = ^u n e n . 

n=l 

If it is known that this function has a branch point, we can try to transform 
the original series (|1 14j) in an implicit function 

F(e,f)=0, 

and determine all required branches of it. 

For this purpose we construct a polynomial F p (e, f) of degree p > 2 

p m 
m=l k=0 

It was assumed Co.i = 1, and the remaining coefficients must be deter- 
mined from the condition 

F p (eJ N (e))=0(e N + 1 ) at e -+ 0. (115) 

Polynomial F p contains 0.5 (p 2 ± 3p — 2) unknowns, the condition (|115[) 
yields N linear algebraic equations, therefore, should N — 0.5 (p 2 ± 3p — 2). 
Once the polynomial F p is found, one can easily find p branches of the 
solution from the equation 

F v = 0. 
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For t he analysis of bifurcations of th ese solutions one can use Newton's 



polygon (|Vainberg and Trenoginl [1974). If a priori information about the 



searching function is known, it can be taken into account for constructing 
the polynomial F p . 

4.5 Estimates of Effective Characteristics of Composite Materi- 
als 

We consider a macroscopically isotropic 2D composite material consist- 
ing of a matrix with inclusions. The aim is to determine the effective con- 
ductivity q from the known matrix (qi) a nd inclusions ((72) conductivitie s 

and the volume fraction (cp). As shown in ( Tokarzewski and Telegal fl997 ), 

£2 

Qi 

conductivity can be written as follows 



if we take e = | - 1 as a small parameter value, the required effective 



— = l + (pe- 0.5(^(1 - cp)e 2 + 0(e 3 ). (116) 



Using the first two terms of expansion (|116[L one obtains 

< — < 1 + tpe, 



l + <pe qi 
hence the bounds of Wiener 
1-tp (p 



Qi Q2 



1 

< q < (1 - cp)q 1 + (pq 2 . 



4.6 Continualization 



We study a chain of n + 2 material points with the same masses m, 
located in equilibrium states in the points of the axis x with coordinates 
jh(j = 0, 1 . . . . , n, n + 1) and susp ended by elastic couplings of stiffness c 



jn{j = u, 1 . . . . , n, n + i J ana susp 
(Figure 6) (jAndrianov et all l2Q10h 




Figure 6. A chain of elastically coupled masses. 
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Owing to the Hooke's law the elastic force acting on the j-th mass is as 
follows 

0jCO = c[y j+1 (t) - yj(t)] - c[yj(t) - yj-i(t)] = c[ % _i(t) - 2y j (t) + y j+ i(t)], 

where j = 1, 2, . . . , n and yj (t) is the displacement of the j-th material point 
from its static equilibrium position. 

Applying the second Newton's law one obtains the following system of 
ODEs governing chain dynamics 

ma jtt (t) =c{(T j+1 -2(jj +aj-i), j = 1,2, ...,n. (117) 

Let us suppose the following BCs 

<ro(t) = a n+1 (t)=0. (118) 

For large values of n usually continuum approximation of discrete problem 
is applied. In our case it takes the form of 

ma u (x,t) = ch 2 cr xx (x,t), (119) 
(j(0, t) = a(e,t) = Q. (120) 

Formally, one can rewrite Eq. (|117p as pseudo-differential equation: 

d 2 a A . 2 / ih d \ „ 
™— + 4csm^- y -j. = 0. (I2l) 

Pseudo-differential operator can be split into the Maclaurin series as follows 

o / ih d\ l „ h 2k d 2k 

sm — T«Z = -o E 



2 fey 2 fc f i (2fc)! &r 2fc 



ft 2 a 2 / h 2 d 2 h A d A h 6 d 6 
l + — - 



4 dx 2 V 12 &r 2 360 fe 4 10080 dx 6 

(122) 

With only keeping the first term in the last line of Eq. (|122p , one obtains 
a continuous approximation (|119p . Keeping the first three terms in Eq. 
(|122p . the following model is obtained 

d 2 a u2 ( d 2 h 2 d A h A d 6 \ 
m W =Ch (^ + 12^ + 360^J (123) 

In the case of periodic BCs for a discrete chain one obtains the following 
BCs for Eq. fiTISh : 

cr = (J X x = (?xxxx = for x = 0, L (124) 
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BVP (j!23p , (|124j) is of the 6th order in spatial variable. Using PA we can 
obtain a modified continuous approximation of the 2nd order. If only two 
terms are left in in the last line of Eq. (|122p . then the PA can be cast into 
the following form: 

d 2 h 2 d 4 _ dx 2 
dx 2 12 dx 4 ~ h 2 d 2 ' 
~ 12 &^ 

For justification of this procedure Fourier or Laplace transforms can be 
used. 

The corresponding so-called quasicontinuum model reads 

mil- j^-q^2) <*tt - ch 2 a xx = (125) 
The BCs for Eq. (fT25l) have the form (fT20j) . 

4.7 Rational interpolation 

Here we follow ([Floater and Hormanni 2007 ). A simple way to approxi- 



mate a function is to choose a sequence of points 

a = xo < x\ < x 2 • • • < x n = 6, 
and to construct the interpolating polynomial p n (x) 

Pn(xi) = f(xi), i = 0,1,2,..., n. 

However, as is well-known p n (x) may not be a good approximation to /, 
and for large n ^> 1 it can exhibit wild oscillations. If we are free to choose 
the distribution of the interpolation points x^, one remedy is to cluster them 
near the end-points of the interval [a, 6], for example using various kinds of 
Chebyshev points. 

Sowe have to make do with them, and then we need to look for other 
kinds of interpolants. A very popular alternative nowadays is to use splines 
(piecewise polynomials), which have become a standard tool for many kinds 
of interpolation and approximation algorithms, and for geometric modeling. 
However, it has been known for a long time that the use of rational functions 
can also lead to much better approximations than ordinary polynomials. In 
fact, both, polynomial and rational interpolation, can exhibit exponential 
convergence when approximating analytic functions. In "classical" rational 
interpolation, one chooses some M and TV such that M + N = n and fits 
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a rational function of the form to the values where pm and 

are polynomials of degrees M and N respectively. If n is even, it is typical 
to set M + N — ^, and some authors have reported excellent results. The 
main drawback, though, is that there is no control over the occurrence of 
po les in the interval of inte r polati on. 

Berrut and Mittelmannl ( 1997 ) suggested that it might be possible to 



avoid poles by using rational functions of higher degree. They considered 
algorithms which fit rational functions whose numerator and denominator 
degrees can both be as high as n. This is a convenient class of rational inter- 
polants because every such interpolant can be written in so-called barycen- 
tric form 

n \. 

E -*-f{xi) 

-0 X — Xi 



r(x) — 



n \ 

E -z-ffr 

i=0 X Xi 



for some real values A^. Thus it suffices to choose the weights A^ in order 
to specify r, and the idea is to search for weights which give interpolant s r 
that have no poles and p referably good approximation properties. Various 



approach is described in iFloater and Hormannl (|2007l ), in particular, one 



can choose A^ = (—1)% i = 0, 1, 2, . . . , n. 

4.8 Some Other Applications 

PA is widely used for the construction of solitons and other localized 
solutions o f nonlinear problems, even in conne ction with the appeared term 
"padeon" ([Lambert and Musette , I1984L Il986h , As a simple model, we con- 



sider the nonlinear BVP 

y"-y + 2y 3 = 0, (126) 
2/(0) = 1, y(oo)=0, (127) 

that has an exact localized solution 

y = cosh" 1 0) (128) 

Quasilinear asymptotics give a solution in the following form: 

y = Ce~ x (1 - 0.2bC 2 e~ 2x + 0.0625C 4 e - 4 * + ...), C = const. (129) 

It is easy to verify that with reconstructing the truncated series (j!29j) in the 
PA and to determine the constant C from the BCs (|127|) and we arrive at 
the exact solution (|128|) . 
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It is also interesting to use the PA to the problems with the phenomenon 
of "blow-up", when the solution goes to infinity at a finite value of the 
argument. For example, the Cauchy problem 



dx o 
— = ax + ex , 
at 



x(Q) = 1, 0<£«a«l, 



(130) 



has the exact solution 



x(t) = 



aex.p(at) 



a + e — eex.p(at) 



(131) 



which tends to infinity for t — »• In [(a + s)/e\. 
Regular asymptotic expansion 

x(t) ~ exp(at) — ea~ x exp(a£)[l — exp(at)] + . . . 

can not describe this phenomenon, but the use of the PA gives the exact 
solution (IT3T1) . 

PA allows to expand the scope of the known approximate methods. For 
example, in the method of harmonic balance the representation of the solu- 
tion of a rational function of the type 



x(t) 



N 

{A n cos[(2n + l)ut] + B n sin[(2n + l)ut]} 

n=0 

N 

1+5^ {Cm cos(2mut) + Dm s'm(2mujt)} 

n=0 



substantially increases the accuracy of approximation (Handy, 1985; Mickens, 
1986). PA can be used effectively to solve ill-posed problem s. This could in- 
clude recon struction of functions in the presence of noise (jGilewicz and Pindor , 
1997L Il999h , various problems of dehomogenization (i.e., determining the 



components of a comp osite material on its homogenized c haracteristics) 
(ICherkaev and Ou. ,2008), etc. We must also mention 2D PA (Vavilov et al 



2002). 



5 Matching of Limiting Asymptotic Expansions 

5.1 Method of Asymptotically Equivalent Functions for Inver- 
sion of Laplace Transform 

This method was originally proposed by Slepian and Yakovlev for the 
treatment of integral transformations. Here is a description of this method, 
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following I Slepvan and Yakovlevl (|l980l ). Suppose that the Laplace transform 
of a function of a real variable f(t) is: 



oo 

F(s) = J f(t)e~ st ds. 



To obtain an approximate expression for the inverse transform, it is 
necessary to clarify the behavior of the transform the vicinity of the points 
5 = and s = oo and determine the nature and location of its singular 
points are on the exact boundary of the regularity or near it. Then the 
transform F(s) replaced by the function Fq(s), allowing the exact inversion 
and satisfying the following conditions: 

1. Functions Fo(s) and F(s) are asymptotically equivalent at s — » oo and 
s — » 0, i.e. 

Fo(s) ~ F(s) at s — ^ and s — >> oo. 

2. Singular points of the functions Fo(s) and F(s), located on the exact 
boundary of the regularity, coincide. 

The free parameters of the function Fq(s) are chosen in such a way that 
they satisfy the conditions of the approximate approximation of F(s) in the 
sense of minimum relative error for all real values s > 0: 



mm < max 



Fo(s,oti,ot2, • • • ,otk) _ 1 



F(s) 



(132) 



Condition (|132p is achieved by variation of free the parameters ol{. Often 
the implementation of equalities 



oo oo 

J F (s) ds = J F(s) ds 



or Fq ~ Fq at s —> leads to a rather precise fulfillment of the requirements 
(|132p . Constructed in such a way the function is called asymptotically 
equivalent function (AEF). 

Here is an example of constructing AEF. Find t he inverse transform if the 



Lapla ce transform is the modified Bessel function (|Abramowitz and Stegun . 



19651 . Sect. 9): 

c 

K (s) = - ln( S /2)7 ( S ) + £ ^rr^*(k + 1) (133) 



g 2k 



2 2k (k\) 2 



40 



where ^(z) is the psi (digamma) function (jAbramowitz and SteeurJ . Il965 . 
Sect. 6). 

For pure imaginary values of the argument s(s = iy; < \y\ < oo) 
function Kq(s) has no singular points. Consequently, we can restrict the 
study of its beh avior for s — > and s — > oo. Th e corresponding asymptotic 
expressions are ( Abramowitz and Steeunl . Il965. Sect. 9): 



K (a) 
K (s) 



ln-+ 7 _ 




O(s), 



s -)■ 0, 



(134) 



OO, 



where 7 is the Euler's constant (7 = 1 , 781 . . . ) (note the typo in the first 
formula (fT3l) in lSlepvan and Yakovlevl (|l980h ). 

The analyzed Laplace transform has a branch point of the logarithmic 
type, branch point of an algebraic type, and an essential singularity. These 
singular points need to be stored in the structure of the zero approximation. 
The most simple way to obtain such a structure, combining two asymptotic 
representations (|134p so that they are mutually do not distort each other 
and contain free parameters, which could be disposed of in the future. As 
a result, we arrive at the zero approximation 



Fo(s) 



ln- 



1 



2 y/ITP. 



(135) 



where a and j3 are the free parameters. 

It is easy to see that expression (|135p has the correct asymptotic be- 
havior s — »> 00. The free parameters are determined from the condition of 
coincidence of the asymptotics of the functions Ko(s) and Fo(s) for s — >• 
and the equality of integrals 



00 00 
J F (s) ds = J K (s) ds. 



As a result of calculations one obtains a system of transcendental equa- 
tions 



In a 



ln2- 7 ; 



In a — e a Ei(— a) + 7 ■ 



[l-erf(V2)] 



7T 

2' 
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where Ei(. . . ) is the the sine integral (I Abramowitz and Steguni 19651 Sect. 5), 
erf(. . . ) is the the error functi on (Abramowitz and SteffunTll965l Sect. 7) 
(note typo in these formulas in iSlepvan and Yakovlevl (|198Q[ )). 

Solving the prescription system numerically, one finds a = 0.3192, f3 = 
0.9927. 

Then the approximate inverse transform can be written as follows: 



fo(t) 



1 - exp[-a(i - 1)] exp[-0(t - 1)] 



The exact expression for the function f(t) is: 

1 



/(*) = 



V* 5 



:H(t-l) 



H(t-1). (136) 



(137) 



Comparison of exact (|137p (solid line) and approximate (|136|) (dotted line 
with circles) inversions is shown in Figure 7. As it can be seen, a satisfactory 
result is obtained even in the zero approximation. Analogously, one can 



f(t) 




0123456789 10 

Figure 7. Comparison of the exact Laplace transform inversion with the 
treatment by the method of AEFs. 



construct AEFs for inverse sine and cosine Fourier transforms, Hankel and 
other integral transforms. 
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5.2 Two-point Pade Approximants 

The analysis of numerous examples confirms: usually implemented a sort 
of "complementarity principle": if for e —> one can construct a physically 
meaningful asymptotics, there is a nontrivial asymptotics and e — >• oo. The 
most difficult in terms of the asymptotic approach is the intermediate case 
of e ~ 1. In this domain numerical methods typically work well, however, if 
the task is to investigate the solution depending on the parameter e, then it 
is inconvenient to use different solutions in different areas. Construction of 
a unified solution on the basis of limiting asymptotics is not a trivial task, 
which can be summarized as follows: we know the behavior of functions in 
zones I and III (Figure 8), we need to construct it in the zone II. For this 
purpose one can us e a two-point Pade approximant s (TPPA). We give the 
definition following iBaker and Graves- Morris! ( 19961 ). Let 



F(e) = °i £i at £ °> ( 138 ) 

2=0 

OO 

F(e) = c i e ~ l at e -> 00 ■ ( 139 ) 



i=0 
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Its TPPA is a rational function of the form 



a + aie ■ 

f[n/m}{£) = YTh~e 



■ a n e" 



!£ + ••• + b m e m 

with k < m + n — 1 coefficients which are determined from the condition 

(1 + bis H h bme 171 ) (co + cie + c 2 £ 2 + . . . ) = a + aie H h a n e n 

and the remaining m + n — k coefficients of a similar condition for e -1 . 

As an example of TPPA using for matching of limiting asymptotics, 
consider the solution of the van der Pol equation: 

x + ex (x 2 — l) + x = 0. 

Asympto tic exp r ession s of the oscillation period for small and large val- 
ues of e are (|Hinchl . ll99lh : 



2tt 1 



16 3072 



5£ 4 \ 



at e —> 0, 



T = e(3-21n2) at e 



oo. 



(140) 



(141) 



For constructing TPPA we use the four conditions at e — » and the two 
conditions at e —> oo, then 



where 



ao = 27r; ai 



T(e) 



tt 2 (3- 2 In 2) 



a + aie + a 2 e 2 + «3^ 3 
l + M + 6 2 e 2 ' 

tt(3 - 21n2) 2 



(142) 



4(3 - 2 In 2) 2 -7T 2 



a 2 



a 3 

&2 = 



7r 2 (3-21n2) 



16 [4(3-21n2) 2 -7r 2 ] 

tt 2 

16 [4(3-21n2) 2 -7r 2 ]' 



2 [4(3 -2 In 2) 2 - tt 2 ] 
7r(3-21n2) 
2 [4(3 -2 In 2) 2 - tt 2 ] ' 



Table 4 show s the results of the com parison of numerical values of the pe- 
riod, given in I Andersen et al.l (|1984f ), with the results calculated by formula 
(fH2l) . 

Now we construct inverse Laplace transform with the TPPA. Let the 
original function is as follows: 



/(t) = (i + 



2\-0.5 



(143) 
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£ 


T numerical 


T numerical 


1 


6.66 


6.61 


2 


7.63 


7.37 


3 


8.86 


8.40 


4 


10.20 


9.55 


5 


11.61 


10.81 


6 


13.06 


12.15 


7 


14.54 


13.54 


8 


16.04 


14.96 


9 


17.55 


16.42 


10 


19.08 


17.89 


20 


34.68 


33.30 


30 


50.54 


49.13 


40 


66.50 


65.10 


50 


82.51 


81.14 


60 


98.54 


97.20 


70 


114.60 


113.29 


80 


130.67 


129.40 


90 


146.75 


145.49 


100 


162.84 


161.61 



Table 4. Comparison of numerical results and calculations using the TPPA. 



Asymptotics of this function looks like: 



/(*) 



1 - 0.5t 2 
t' 1 + ... 



at t -> 0, 
at t —> oo. 



TPPA in this case can be written as: 

/w i+a5f 



l + 0.5i + 0.5i 2 ' 



(144) 



Numerical results are shown in Figure 9. An approximate inversion (| 144ft 
(upper curve) agrees well with the original (|143p (lower curve) for all values 
of the argument. 



Oth er example on the effect i ve use of the TPPA s ee |A ndriano v and Awreicewicz 



(200lh; 



19761) : 



Andrianov et al. ( 20041 ): Awreicewicz et al. ( 1998f k lFrost and Harperl 



Wenigeil dl98flh 
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f(t) 

0.8 
0.6 
0.4 
0.2 











































2 A 
Figure 9. Exact 


\ 6 t 8 

and approximate 



5.3 Other Methods of AEFs Constructiong 

Unfortunately, the situations where both asymptotic limits have the form 
of power expansions are rarely encountered in practice, so we have to resort 
to other methods of AEFs constructing. Consider, for example, the BVP 

ey xx -xy = ey, 2/(0) = 1, y(oo) = 0,e «1. (145) 

Solution for small values of x can be written as follows: 

y = l-at+le + 0(e) (146) 



where £ = xe 3 ? a is the arbitrary constant. 

The solutio n for large values of x is constructed using the WKB method 
(jNavfehl . l2000h 



y = b£ 4 exp 
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(147) 



where b is an arbitrary constant. 

Now we match these asymptotics. Because of the exponential in Eq. 
(|147p using TPPA in the original form is not possible. Therefore, we con- 
struct AEF, based on the following considerations: for large values of the 
variable £ exponent from Eq. (|147|) is taken into account in its original 
form, and for small values of the variable £ it is expanded in a Maclaurin 
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series. Constructed in this way AEF has the form 

2^3 2^5 32 



4 



Va = ; 32 a ^ eXp ' (M8) 

1 + y^ 4 



2 



The coefficients a and b in Eq. (|148j) still remain uncertain. For calculation 
of these constants one can use some integral relations, for example, obtained 
from Eq. (|148|) by multiplying them with the weighting functions 1, x, x 2 , . . . 
and further integration over the interval [0,oo). In the end, such values of 
the constants are found: 



3/ , V3 



2^ ' 



(149) 



Numerical calculations show that the formula (|148[) with constants (|149p 
approximates the desired solution in the whole interval [0, oo) with an error 
not exceeding 1.5%. 

When choosing the constants one can use other methods, then a lot 
depends on the skill of the researcher. Of course, it is necessary to ensure 
the correct qualitative behavior of AEFs, avoiding, for example, do not 
correspond to the problem of zeros of the denominator. To do this, one 
can vary the number of terms in the asymptotics and the numerator and 
denominator constructed uniformly suitable solutions. In general form the 



method of rational AEF can be described as follows (Martin and Baker 



1991f ). Let us assume that function f(z) has the following asymptotics: 

f(z) = F(z) at z -> oo, (150) 

and 

CO 

f(z) = c i zi at z 0- (151) 

i=0 

Then the AEF can be produced from the Eqs. (|152p . (|153p as follows: 

m 

/(*)»^ at z^O. (152) 

E 

where a^, Pi are considered not as constants but as some functions of z. 
Functions ai{z) and f3%{z) are chosen in such a way that: 
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1. the expansion of AEF (j!52p in powers of z for z — >• matches the 
perturbation expansion (|151|) : 

2. the asymptotic behavior of AEF (|152p for z — >• oo coincides with the 
function F(z) (fT50|) . 

In the construction of AEFs a priori qualitative information is very im- 
portant. For example, if from any considerations it is known that the un- 
known function is c lose to the power, you can use the method of Sommerfeld 
(jKashin et al.Ul989h . Its essence is to replace a segment of the power series 



to the function 



f(x) = 1 + a x x + a 2 x z + • • • , (153) 



f{x) « (1 + Axy. (154) 



Expanding expression (|154|) in a Maclaurin series and comparing coefficients 
of this expansion with (|153p , one obtains 

a\ — 2d2 a\ 



ai a\ — 2a2 

Numerical approaches also can be used for the construction of AEFs. 
For example, in paper by iKaas- Petersen! (|l987[ ) a computational technique 



for matching limiting asymptotics is described. 

Sometimes it is possible to construct so called composite equations, 
which can be treated as "asymptotically equival ent equati o ns". L et us em- 
phasize, that the composite equations, due to Ivan Dvkd (|l975bh . can be 



obtained in result of synthesis of the limiting cases. The principal idea of 
the method of t he com posite equations can be formulated in the following 
way (jvan Dvkd . fl975bl p.195): 



1. Identify the terms in the differential equations whose neglect in the 
straightforward approximation is responsible for the nonuniformity. 

2. Approximate those terms insofar as possible while retaining their es- 
sential character in the region of nonuniformity. 

Let's dwell on the terminology. Here we use the term "asymptotically 
equivalent fu nction". Other term s ("reduced method of matched a symptotic 



expansions" (Kash in et a" 



,, 1989 ), "quasifr actional app roximants" ( Chalbaud and Martini . 



expj_ 

Il992l ). "mimic function" ([Gaunt and Guttmanl ll974L p.181-243) also used. 



5.4 Example: Schrodinger Equation 

For the Schrodinger equation ([60]) with boundary conditions (|6Tj) previ- 
ously we obtained a solution for the exponent, little different from the two 
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([68]). In ( Boettcher and Benderl Il990l ) the following asymptotic solutions 
for TV — » oo is obtained: 



E (N) = ^(2N)-^T (JL. 



(155) 



Using (|68j) and (|155p , we construct AEF 



7r + r 



^o(AT) 



TV 



TV + 1 



4(27V + a)^H 



(156) 



where a = 7r 2 r(1.25) — 2 « 6.946. Numerical results are presented in Table 
5. It is evident that formula (jl56j) gives good results for all the values of N. 



N 


Eq numerical 
(Boettcher and Bender, 
(1990) 


Ea.<nm> 


Error,% 


1 


1.0000 


1.0 





2 


1.0604 


0.9974 


5.9364 


4 


1.2258 


1.17446 


4.1882 


10 


1.5605 


1.5398 


1.33 


50 


1.1052 


2.1035 


0.079 


200 


2.3379 


2.3376 


0.006 


500 


2.4058 


2.4058 


« 


1500 


2.4431 


2.4431 


« 


3500 


2.4558 


2.45558 


« 



Table 5. Comparison of numerical and analytical results of the energy levels 
for the SchrAudinger equation. 



5.5 Example: AEFs in the Theory of Composites 

Now let us consider an application of the method of AEFs for the cal- 
culation of the effective heat conductivity of an infinite regular array of 
perfectly conducting s phere s, embedded in a matrix with unit conductivity. 
Sangani and Acrivosl (p. 983) have obtained the following expansion for the 



49 



effective conductivity (k): 
(k) 



1 



3c 



-1 



a\c 3 



1 + a 2 c s 

a 3 d 



a^c 3 -f a^c 3 



a^c 3 



O 



(of) 



(157) 

where c is the volume fracture of inclusions. Here we consider three types of 
space arrangement of spheres, namely, the simple cubic (SC), body centered 
cubic (BCC) and face centered cubic (FCC) arrays. The constants a$ for 
these arrays are given in Table 6. 





CLi 


«2 


«3 


CL4 


a 5 


a 6 


SC array 


1.305 


0.231 


0.405 


0.0723 


0.153 


0.0105 


BCC array 


0.129 


-0.413 


0.764 


0.257 


0.0113 


0.00562 


FCC array 


0.0753 


0.697 


-07.41 


0.0420 


0.0231 


9.14- 10- 7 



^max 



Table 6. The constants ai, . . . , a§ in Eq. (|157|) . 

In the case of perfectly conducting large spheres (c — » c maxi where c„ 
is the maximum volume fraction for a sphere) the problem can be solved 
by means of a reasonable physical assumption that the heat flux occurs 
entirely in the region where spheres are in a near contact. Thus, the effective 
conductivity is determined in the asymptotic form for the flux between two 
spheres, which is logarithmically singular in th e width of a gap, justifying 
the assumption ([McPhedran and Milton , 198lh : 



(k) = -M 1 ln X -M 2 + 0( X ~ 1 ), (158) 
where y = 1 — ( — — ) 3 is the dimensionless width of a gap between the 

y Cmax J 

neighboring spheres, x ~~ ^ f° r c ~~ ^ c m, Mi = 0.5c max p, p is the number of 
contact points at the surface of a sphere; M2 is a constant, depending on 
the type of space arrangement of spheres. The values of Mi, M2 and c max 
for the three types of cubic arrays are given in Table 7. 

On the basis of limiting solutions (|157[) and (|158[) we develop the AEF 
valid for all values of the volume fraction of inclusions c G [0,c ma J: 

_ flM + jW-P + flh, (159) 
Q{c) 

Here the functions Pi(c), Q(c) and the constants P2, P3 are determined as 
follows: 
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Mi 


M 2 




SC array 


tt/2 


0.7 


tt/6 


BCC array 


VStt/2 


2.4 




FCC array 




7.1 





12 



Table 7. The constants Mi, M 2 and 




0.5 C 



Figure 10. Effective conductivity (k) /k m of the SC array vs. volume 
fraction of inclusions c. 



<b/k n 




0.1 0.2 0.3 0.4 0.5 0.6 C 



Figure 11. Effective conductivity (k) /k m of the BCC array vs. volume 
fraction of inclusions c. 
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<hik m 




0.10.2 0.3 0.4 0.50.60.7c 



Figure 12. Effective conductivity (k) /k m of the FCC array vs. volume 
fraction of inclusions c. 



Q(c) = 1 — c — aiC3° , Pi(c)=^aics, P 2 = for n = 1, 

2=0 

) + Q( 

The AEF (|159[) takes into account leading terms of expansion ([157)) and 
leading terms of expansion (|158[) . Coefficients are: 

1 n Q(c m ax)Mi Q(c ma ^)Mi 

a = 1, « 3 = 2 , a i0 = ai To > 

^C ma:c i n 3 

±\J^max 

Q{ c max)Mi . . , 

otj = J , J = 1,2, ...,ra- Ira, , j 7^ 3, 10. 

The increment of ra and n leads to the growth of the accuracy of the 
obtained solution (j!59j) . Let us illustrate this dependence in the case of 
SC array. We calculated (k) for different values of m and n. In the Figure 
10 o ur analytical results a re co mpared with experimental measurements 
from Meredith and Tobiasl (Il960l) f b l ack d ots). Details of these data can be 
found in McPhedran and McKenziel ( 19781 ). Finally, we restrict ra = 19 and 



n = 2 for all types of arrays, as they provide a satisfactory agreement with 
numerical datas and a rather simple analytical form of the AEF (|159p . 

Numerical results for the BCC and the FCC arrays are displayed in 
Figures 11 and 12 respectively. For BBC array the obtained AEF (|159|) 
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is com pare d with the experimenta l results from iMcKenzie and McPhedran 



( 19771 ) and IMcKenzie et al.l (| 19781 ). For FCC array the experimental data 
are not ava ilable, therefore we are comparing with the numerical results 
obtained by IMcKenzie et aD ( 1978 ) using the Rayleigh method. The agree- 
ment between the analytical solution (|159p and the numerical results is quite 
satisfactory. 
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Institute of General Mechanics, RWTH Aachen University, D-52062 Aachen , 

Germany 

Abstract In this lectures various methods which give a possibility 
to extend an area of applicability of perturbation series and hence to 
omit their local character are analysed. While applying asymptotic 
methods as a rule the following situation appears: the existence of 
asymptotics for e — >• implies an existence of the asymptotics for 
e — >• oo. Therefore, the idea to construct one function valid for the 
whole parameter interval for e is very attractive. The construction 
of asymptotically equivalent functions possessing a known asymp- 
totic behaviour for e — >• and e — >• oo will be discussed. Using 
summation and interpolation procedures we focus on continuous 
models derived from a discrete micro- structure. Various contin- 
ualization procedures that take the non-local interaction between 
variables of the discrete media into account are analysed. 
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1 Overview of Manuscript Preparation and Delivery 

Asymptotic analysis is a constantly growing branch of mathematics which 
influences the development of various pure and applied sciences. The fa- 
mous mathematicians ? and ? said that an asymptotic description is not 
only a suitable instrument for the mathematical analysis of nature but that 
it also has an additional deeper intrinsic meaning, and that the asymptotic 
approach is more than just a mathematical technique; it plays a rather 
fundamental role in science. And here it appears that the many existing 
asymptotic methods comprise a set of approaches that in some way belong 
rather to art than to science. ? even introduced special term "asymptotol- 
ogy" and defined it as the "art of dealing with applied mathematical systems 
in limiting cases". Here it should be noted that he called for a formalization 
of the accumulated experience to convert the art of asymptotology into a 
science of asymptotology. 

Asymptotic methods for solving mechanical and physical problems have 
been developed by many authors. We can mentioned the excellent mono- 
graphs by ?, ?, ?, ??, ??, ?, and many others. The main feature of the 
present book can be formulated as follows: it deals with new trends and 
applications of asymptotic approaches in the fields of Nonlinear Mechan- 
ics and Mechanics of Solids. It illuminates the developments in the field of 
asymptotic mathematics from different viewpoints, reflecting the field's mul- 
tidisciplinary nature. The choice of topics reflects the authors' own research 
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experience and participation in applications. The authors have paid special 
attention to examples and discussions of results, and have tried to avoid 
burying the central ideas in formalism, notations, and technical details. 

2 Some Nonstandard Perturbation Procedures 

2.1 Choice of Small Parameters 

The choice of the asymptotic method and the introduction of small di- 
mensionless parameters in the system is very often the most significant and 
informal part of the analytical study of physical problems. This should 
help the experience and intuition, analysis of the physical nature of the 
problem, experimental and numerical results. It often dictated by physical 
considerations, is clearly shown as dimensionless and scaling. However, it is 
sometimes advantageous to use is not obvious, and perhaps even strange at 
first glance, the initial approximation. To illustrate this, consider a simple 
example (?): algebraic equation 

x 5 +x = l. (1) 

We seek the real root of Eq. (pQ) , the exact value of which can be determined 
numerically: x = 0.75487767... . A small parameter e is not included 
explicitly in Eq. (pQ). Consider the various possibilities of introducing a 
parameter e in Eq. ([1]). 

1. We introduce a small parameter e with a nonlinear term of Eq. (pQ) 

ex 5 +x = l, (2) 

and present x as a series of e 

x = a + a\s + a 2 e 2 + (3) 

Substituting the series (|3]) into Eq. (j2j) and equating terms of equal 
powers, we obtain 

do = 1, a i = —1, a 2 = 5, a 3 = —35, a± = 285, a 5 = —2530, 
a 6 = 23751. 

These values can be obtained by a closed expression for the coeffi- 
cients a n : 

(-l) n (5n)! 
an ~ n!(4n + l)f 

/i 4 

The radius R of convergence of the series © is R = §5 = 0.08192. 
Consequently, for e = 1 series (|3]) diverges very fast, so the sum of 
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the first six terms is 21476. The situation can be corrected by the 
method of Pade approximants (see Sect. 2). Constructing a Pade 
approximant with three terms in the numerator and denominator and 
calculating it with e = 1, we obtain the value of the root x — 0.76369 
(the deviation from the exact value is 1.2%). 
2. We now introduce a small parameter e near the linear term of Eq. (pQ) 

£ 5 +££ = l. (4) 

Presenting the solution of Eq. (|4]) in the form 

x(e) = b + he + b 2 e 2 + . . . , (5) 
we have, after applying the standard procedure of perturbation method 



b = 1, h = -1, b 2 = -I05, cl 3 = -Itof, b A = 0, b 5 = 
21 h 78 Ib lIb 

15625' U6 ~ 78125* 

And in this case we can construct a general expression for the co- 
efficients with: 

r[(4n-l)/5] 
5r[(4-n)/5]n!' 

and determine the radius of convergence of the series R = 
4(4 5 /5) = 1, 64938 .... The value of x(l), taking into account the first 
six terms of the series (j5j), deviates from the exact by 0.07%. 
3. We introduce the "small parameter" 5 in the exponent 

x 1+5 +x = l, (6) 

and represent x in the form 

x = c + c\ 5 + c 2 S 2 + (7) 

In addition, we use the expansion: 

x 1+5 =x(l + (51n|x| + ...). 

Coefficients of the series ([7]) are determined easily: 

c = 0.5, ci=0.251n2, c 2 = -0.125 In 2, .... 

The radius of convergence is one in this case, and calculate when it 
should be 5 = 4. Using Pade approximants with three terms in the nu- 
merator and denominator, if e — 1, we find x = 0.75448 that only devi- 
ates from the exact by 0.05% result. Calculating q for i = 0, 1, . . . , 12 
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and constructing Pade approximant with six terms in the numera- 
tor and denominator, we find x = 0.75487654 (0.00015% error). The 
method is called "the method of small delta" (see Sect. 12. 3p . 
4. We now assume the exponent as a large parameter. Consider the 
equation 

x n + x = l. (8) 



Assuming n —> oo (a method of large J, see Sect. 12. 4[L we represent 
the desired solution in the form 



x = 



-(l + zi +z 2 + ...) 



l/rj 



(9) 



where 1 > x\ > X2 > 

Substituting the ansatz (0 in Eq. flS) and taking into account that 



1 



1 



Inn 



1 + - ln(l + x x + x 2 + ...) + ... . 

n 



one obtains in order of increasing accuracy of the formula 

m n x 

n 



In n — In In n 



l/n 



(10) 

(ii) 



For n = 2 formula (pTUj) gives x = 0.58871; the error compared to the 
exact solution (0.5(v / 5 - 1) « 0.618034) is 4.7%. When n = 5, we 
obtain x = 0.79715 from Eq. ([TU|) (from the numerical solution one 
obtains x = 75488; error 4.4%). For n = 5, Eq. (11]) gives x = 0.74318 
(error 1.5%). Thus, even the first terms of the large S asymptotics give 
excellent results. 

Hence, in this case the method of large delta provides already a 
good accuracy at low orders of perturbation method. Approximation 
([T0|h ([TT]) give an example of nonpower asymptotics. 



2.2 Homotopy Perturbation Method 

In recent years the so-called homotopy perturbation method was pop- 
ular (?) (term "method of artificial small parameters" is also used). Its 
essence is as follows. In the equations or boundary conditions a parameter 
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s is introduced so that for e = one obtains a BVP which admits a simple 
solution, and for e = 1 one obtains the governing BVP. Then the pertur- 
bation method to e is apllied and in the resulting solution was adopted. 
Of course, it is not new and was already used by ?, chapt. 1.6 ? and ?. 
The novelty (and, in my opinion, successful) is the title, emphasizing the 
continuous transition from the initial value e = to the value of e = 1 (ho- 
motopy deformation). Let us analyse an example of homotopy perturbation 
parameter method using ?. A special feature of nonlinear systems with dis- 
tributed parameters is the possibility of internal resonance between modes. 
That is why in many cases the neglection of higher modes can lead to sig- 
nificant errors. The following describes the asymptotic method of solving 
problems of nonlinear vibrations of systems with distributed parameters, 
allowing approximately to take into account all modes. The oscillations of 
a square membrane on a nonlinear elastic foundation can be written as: 

d 2 w d 2 w d 2 w o , . 

o^ + W~^- cw - £W °' (12) 

where e is the dimensionless small parameter (e <C 1). 
The BCs are as follows: 

w \ x =o,l = w \y=o,L = °- ( 13 ) 
The desired periodic solution must satisfy the periodicity conditions 

w(t) =w(t + T), (14) 

where T = ^ is the period, and Q is the natural frequency of oscillation. We 
seek the natural frequencies corresponding to these forms of fundamental 
vibration frequencies at which the linear case (e = 0) is realized by one 
half- wave in each direction x and y. We introduce the transformation of 
time: 

r = ojt. (15) 

The solution is sought in the form of expansions 

w = wo + ewi + e 2 w 2 + . . . , (16) 

uj = ujq + euoi + e 2 (jj 2 + (17) 

Substituting ansatzes ([T6]h ([T7|) into Eqs. ([T2]) - jT2J) and equating terms of 
equal powers , we obtain a recurrent sequence of linear BVPs: 

^ + ^--o^--o = 0, (18) 

d 2 w 1 d 2 w 1 2 d 2 w 1 d 2 w 3 , . 

^ + = ^i-^+t^o, (19) 
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The BCs ([T3]) and periodicity conditions ([T4]) take the form for i = 
1,2,...: 

™4=o,l = ^Uo,l = 0, (20) 
^i(r) = iy*(r + 2tt). (21) 
The solution of Eq. (fT5j) is as follows: 



oo oo 



,0 = ^ m ' n sin ( —^—t ) sin (^^) sin (~^y) ' ( 22 ) 

m=l m=l ^ ' 



where cj m?n = y tt 2 ( m2 + n2 ) -f C; m , n = 1, 2, 3,... and A^i is the 
amplitude of the fundamental tone of vibrations; A m?n , m, n = 1, 2, 3..., 
(m, n) 7^ (1,1) is the amplitude of the subsequent modes; 0Jm,n are the 
natural frequencies of the linear system, ujo = uji^. 

The next approximation is a result of solving the BVP ([T9]) - ([2T]) . To pre- 
vent the appearance of secular terms in the r.h.s. of Eq. ([T9]) the coefficients 
of the terms of the form 

/ TTTYl \ / TTT). \ 

1, 2, 3,... 



• ( Um,n \ • (^m \ . [Tin \ 

sin I —r J sin y — j^ x ) sm \~J^y) ' m ' n 



should be equated to zero. 

These conditions lead to an infinite system of nonlinear algebraic equa- 
tions: 

~ , oooooooooooo 

-^f 1 (« m ,») a = EEEEEEC PS) ^4A, S , (23) 

P2 i=l j=l k=l 1=1 p=l s=l 

where m, n = 1, 2, 3, . . . . 

Coefficients are found by substituting ansatz ([22]) into the r.h.s. of Eqs. 
(fT9|) and the relevant simplifications. System ([23]) can be solved by reduc- 
tion. However, a sufficiently large number of retention equations beging to 
show significant computational difficulties. In addition, this approach does 
not take into account the influence of higher oscillation modes. Therefore, 
further we use the homotopy perturbation method. 

On the right side of each (m, n)-th equation of system ([23]) we introduce 
the parameter fi to those members of AijAkjA p ^ s , for which the following 
condition is valid: (i > m) U (k > m) U (p > m) U (j > n) U (I > n) U (s > n). 
Thus, for \i = system ([23]) takes the "triangular" appearance, and for 
/i = 1 it returns to its original form. Next, we seek a solution in the form 
of expansions: 

wi = 10^+^+^10^+..., (24) 

A m ,n = ^-m]n + f lj ^m,n + f 1 ' -^m.n + • • • > (25) 
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where ra, n = 1, 2, 3, . . . , (m, n) 7^ (1, 1). 

We put \i — 1 . This approach allows us to keep any number of equations 
in system ([23]) . Below we limit ourselves the first two terms in the expansions 
(O, ([25]) . We analyze the solutions and note that in this problem the 
parameter c plays the role of the bifurcation parameter. In general, for 
c^O, c~ 1 , the system ([23]) admits the following solution: 



hj = 1,2,3,..., ^ (m,n), 

27 A 2 u 
Ul Too ' ?™,n= 1,2,3, .... 

Amplitude-frequency response is given by: 

A 2 

On,n = Um,n + 0.2109375^^ £ + . . . . 

Of particular interest is the case when the linear component of the restor- 
ing force is zero (c = 0) and the phenomenon of internal resonance between 
the modes of oscillations occurs. Solving system ([23]) by the described 
method, we find 

A m ,n = 0, ra, n = 1,2,3,..., (m, n) ^ (1, 1), (m, n) ^ (2i - 1, 2i - 1), 
i = 1, 2, 3 . . . , A 3 , 3 = -4.5662 • l(r 3 A M , A 5 , 5 = 2.1139 • 1(T 5 A M , . . . , 
cji = 0.211048A? )1 /o;o. 

If the oscillations are excited by the mode (1, 1) all odd modes (3, 3), (5, 5) 
etc. are also realized. However, if the oscillations are excited by one of the 
higher modes, the result of redistribution of energy modes appear at lower 
orders until the fundamental tone (1,1). 

2.3 Method of Small Delta 

? proposed an effective method of small (5, which we explain by examples. 
Let us construct a periodic solution of the following Cauchy problem 

x tt + x 3 = 0, (26) 

x(0) = l, x t (0) = 0. (27) 
We introduce a homotopy parameter S in Eq. ([26]). thus: 

x tt + x 1+25 = 0. (28) 

At the final expression one should put 5 = 1, but in the process of solving 
one must suppose 5 <C 1. Then 

x 26 = 1 + S\nx 2 + O.bS 2 (\nx 2 ) 2 + . . . . (29) 
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We suppose the solution of Eq. (|26|) in the form 

oo 

Y,8 k x k , (30) 



X 

/,'-() 



and carry out the change of independent variable 

t=-, (31) 

LO 

where uj 2 = 1 + a\8 + o^ 2 + . . . . 

The constants (i = 1,2,...) are determined in the problem solving 
process. After substituting ansatzes (|29|) - (|3T|) into Eq. (|28|) and splitting 
it with respect to 5, we obtain the following recurrent sequence of Cauchy 
problems: 

XOrr + X = 0, (32) 

z (0) = l, x 0r (0) = 0; (33) 

xi rr + a?i = -x ln(^o) - aixorr, (34) 

xi(0)=xi r =0; (35) 

%2tt + x 2 = —xi hi(x 2 ) - 2xi - x (ln(xo)) 2 - a 2 XQ TT - ol\X\ tt \ (36) 

*i(0)=zi T = 0; (37) 



A Cauchy problem of the zero approximation ([32]) , (|33|) has the solution . 

xo — cos r. 

In a first approximation, one obtains 

x lTT + xi = — cos r ln(cos 2 r) + ql\ cos r = Lo- 

The condition of the absence of the secular terms in the solution of this 
equation can be written as follows: 

tt/2 

J Lo cost dt = 0, 



which is determined by the constant a± = 1 — 2 In 2. The expression for the 
oscillation period can be written as 

T = 2tt[1 + J(ln2 — 0.5)]. 
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For 5 = 1, we have T = 6.8070, while the exact value is T = 7.4164 
(the error of approximate solutions is 8.2%). The solution of the Cauchy 
problems of the next approximation ([36]) , (|37|) gives the value of the period, 
it practically coincides with the exact one (T = 7.4111). 

We now consider the wave equation 

uu = u xx , (38) 

with nonlinear BCs: 

u(0,i) = 0, (39) 

u x (l,t) + u(l,t) + u s (l,t) = 0. (40) 

We introduce the parameter S in Eq. ([40]) as follows: 

1^(1, t) + u(l,t) + u l+25 (l, t) = 0. (41) 

In the final expression should we put 5 = 1, but in the process of decision 
we assume S <C 1. Then 



U 3 EE U 1+25 



5 2 2 

1 + S\nu 2 H (lnw 2 ) 



2 

Suppose further the solution of Eq. ([38]) in the form 



(42) 



J2 sku k- (43) 



IX = 

k=0 

After substituting ansatzes ([43]), (@TJ) into Eqs. ([38]) . ([39l) . (I4T1) and the 
splitting of the parameter 5 we obtain the following recurrent sequence of 
BVPs: 

uorr = u 0xx ; (44) 

at x = 0, u = 0; (45) 

at x = 1, 'Uoa? + 2^o = 0; (46) 
l 

^Orr = Uq xx ^ ^ ^i—pUprr'i (^7) 

at x = 0, ?ii = 0; (48) 
at x = 1, i^i^ + 2^i = — uq In i^; (49) 

2 

^Orr = - ^2 a i-P U Prr] (50) 

p=0 
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at x = 0, U2 = 0; (51) 

at x — 1, U2 X + 2^2 = — u\ lniiQ — 2u\ — 0.5^0 ( m ^o) 2 5 (52) 

. . . , 

where ao = 0. 

The solution of the BVP of the zero approximation can be 

written as 

= As'm(uJox) sin(o;oT), (53) 
where the frequency uoq is determined from the transcendental equation 

uo = 2 tanc^o- (54) 
The first few nonzero values of uo are given in Table [IJ 



n (l) 










2.289 


5.087 


8.096 


11.173 


14.276 












,.(10) 


17.393 


20.518 


23.646 


26.778 


29.912 



Table 1. The first few roots of the transcendental Eq. 

When A: — > oo we have the asymptotics: uj^ — 0.57r(2/c + 1). 
The BVP problem of the first approximation are as follows: 

ui xx - ui rr = (XiAuJq sm(uj x) sui(cj t); (55) 

at x = 0, ui = 0; (56) 

at x = 1, ui x + 2^i = Ai sin(cJoT) [ln(A 2 sin 2 ujo) + lnsin 2 (cJoT)] ; (57) 

where Ai = — AsincJo- 

The particular solution of Eq. ([55]) satisfying the BC ([56]) has the form 

= — 0.5a iAujox cos (c^o^) sin (ujqt). (58) 

We choose the constant a\ in such a way that it compensates the secular 
term on the r.h.s. of Eq. (|57|) 

_ 2gi 
cj (6 + cj^) 



11 



where Ri = ln(0.25eA 2 sin 2 u ). 

Non-secular harmonics on the r.h.s. of Eq. ([57)) give the solution 

CO ^ 

= AAi T k sm(oj kx) sm(ojokr) — -, (59) 

k=2 

where = l/[fcu;o cos(kuo) + 2sin(fcu;o)]- 

The complete solution of the first approximation has the form 

(1) , (2) 

Assuming S = 1, we obtain the solution of Eqs. (|38|) - ([4Q|) . 
Let us now consider the Schrodinger equation (?) 

* - x 2N * + E9 = 0, (60) 

tf(±oo) =0. (61) 

Here \I/ is the wave function; is the energy,it plays the role of eigen- 
value. It is shown that the eigenvalue problem (|SDJ) , flST} has a discrete 
countable spectrum E n , n = 0,1,2,... (?). For N = 2 the eigenvalue 
problem (|60]h (|6T|) has an exact solution. Now let TV differ little from 2, 

®xx ~ x 2+25 ^ + £tf = 0. (62) 

Now we use the expansion 

x 26 = l + 5\n(x 2 ) + ..., 

and we will search for the eigenfunction \I/ and the eigenvalue E in the form 
of the following series: 

* = + + S 2 ^ 2 + . . . , (63) 

J5 = £ + SE + £ 2 £ 2 + . . . . (64) 

As a result, after the asymptotic splitting we obtain a recursive sequence of 
eigenvalue problems 

^o xx -x 2 y + £ ^o = 0, (65) 
^ixx ~ x 2 ^i + E ^ 1 + ^x^o = x 2 ^ ln(x 2 ), (66) 



|^i|->0 at \x\^oc, i = 1,2,3,.... (67) 
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The solution of the eigenvalue problem ([65]) , ([67]) has the form 
4 n) = 2n + 1, = e-* 2/2 tf n (x), n = 1, 2, 3, . . . , 

where H n (x) is the Struve function (?, Sect. 12). From the eigenvalue prob- 
lem (|SS|), we find 

jr x 2 e-^^ (x)ln(x 2 )dx 

^(n) _ ^oo 

1 V^2 n ra! 
For n = one obtains Hq(x) = 1 and 

oo 

J x 2 lnxe- x2 dx= ^(2-21n2-C), 

— CO 

where C = 0.577215 ... is the Euler constant. Hence 

4 o) = i+4( 2 - 2m2 - c ) (5 +--- • ( 6s ) 



2.4 Method of Large Delta 

An alternative method of small delta is the method of large delta, which 
we demonstrate in an example of a nonlinear equation 

x tt +x n = 0, n = 3,5,7,... . (69) 

This equation can be integrated with the functions Cs and Sn, introduced 
by Liapunov ? (inversions of incomplete beta functions, ?). Note that 
much later the same (up to normalization) function have been proposed 
by Rosenberg, who called them Ateb- functions ?. However, working with 
these objects is inconvenient, and therefore the problem of the approximate 
analytical solution of Eq. (|69|) in elementary functions arises. We construct 
asymptotics of periodic solutions of Eq. ([69]) at n — » oo. Let the initial 
conditions for Eq. ([69]) be 

x(0) = 0, x(0) = 1. (70) 
The first integral of the Cauchy problem ([69]) , ([70]) can be written as: 



fdx\ _ 2x n+1 



— 



V dt J n + 1 



i-TTTT- ( 71 ) 
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The replacement of the x = A A / 2 , A = 2/(n + 1) , and integration gives us 
a solution in the implicit form 



0<£<l 



After replacing £ = sin A # this implicit solution is transformed into an 
expression that contains a small parameter in the exponent of the integrand: 

o<e<7v/2 

X x/2 t = X J sin- 1+A 6>d<9. 

o 

We now consider the integrand separately: 
sin" 1+A 6> = #" 1+A ( - fl-i+A 



-- Aln— - + .. 

.0 sm 

Expanding this function into a Maclaurin series, one obtains 

sin- 1+A = 6~ l+x + + • • • + 0(A). 

The first term of this expression makes the main contribution of this 
expression, so in the first approximation we can write 

\ x ' 2 t*6 x , i.e. NA 1 ^/*. 

In the original variables one obtains 

arwA- A / 2 sin A (A 1 /V/ A ) (72) 
The solution ([72]) should be used on a quarter-period Twith 

Let us analyze the solution ([72]) , ([73]) . At n = 1 one obtains the exact values 
x = sint, T = 2ir , for n — >• oo one obtains T — >• 4. Expanding the r.h.s. of 
the Eq. ([72]) in a series of t and restricting it to the first term, we obtain 
nonsmooth solution (?). We estimate the error of the solution ([73]) . For 
this we use the expression 

tt/2 

A J sin- 1+A d<9 = 0.5A£(0.5A,0.5) = A u (74) 
o 
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where B(. ..,...) is the beta function (?, Sect. 6). The approximate value of 
the integral on the l.h.s. of Eq. ([72]) is calculated as follows: A2 = (tt/2) x 
. Numerical comparison of the values A\ , A2 and error estimate A is given 
in Tabled 



n 


1 


3 


5 




00 


A! 




1.30 


1.20 




1 


A 2 




1.25 


1.16 




1 


A,% 





5 


3 




~ 



Table 2. Comparison of exact and approximate solutions. 

Thus, the first approximation of the asymptotics for n — ^ 00 already 
gives quite acceptable accuracy for practical purposes, even for not very 
large values of n. Note that the expression ([72]) gives an approximation of 
incomplete beta function (?, Sect. 6) from n = 1 (sinus function) to n = 00 
(linear function). 

2.5 Application of Distributions 

Asymptotic methods are based, generally speaking, on the use of Taylor 
series. In this connection the question arises: what to do with the functions 
of the form exp(— s~ 1 x), which cannot be expanded in a Taylor series in 
e — >• 0, using smooth functions (?). The way out lies in the transition to 
the following distribution (?): 

00 

i7(x)exp(-£- 1 x) = ^(-l) n £ n+ V n )(x), (75) 

n=0 

where S(x) is the Dirac delta function, representing the derivative of the 
Heavisyde function H(x) in the theory of distribution; 5( n \x), n = 1, 2, . . . 
are the derivatives of the delta function. 

We show how formal the formula ([75]) can be obtained. Applying the 
Laplace transform to function exp(— e~ 1 x), one obtains: 

00 

/exp(— £~ 1 x) dx — — - — . 
7 ep + 1 



Expanding the r.h.s. of this equation in a Maclaurin series of e, and then 
calculating inverse transform term by term, one obtains expansion ([75]) . 
Thus, we again use the Taylor series, but now in the dual space. 
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Here is another interesting feature of the approach using distributions: 
a singular perturbated problem can be regarded as a regular perturbated 
one (?). Suppose, for example, 

ey' + y = x > 0, y = 1 at .t = 0. 

This is singular perturbated problem: for e = one obtains a smooth 
solution y = 0, which does not satisfy the given initial condition. However, 
one can seek a solution in the form of a nonsmooth function. Namely, 
assuming z(x) — H(x)y(x), one obtains 

ez' = -z + e6(x). (76) 

The solution of Eq. ([76)) is sought in the form of expansions 

oo 

As a result, one obtains 

z = 0, Z! = 5{x), z n+1 = (-l) n S^ n \ n = l,2,.... (77) 

Note that the expressions (|77]) can go to the smooth functions. To do this, 
it is possible to apply the Laplace transform, then the Pade approximants 
in the dual spase and then to calculate inverse Laplace transform. 

We show other application of the asymptotic method using distribution. 
Consider the equation of the membrane, reinforced with fibers of the small, 
but finite width e. The governing PDE is: 

[1 + 2e$ (y)]u xx + u yy = 0, (78) 

CO 

where $ (y) = ]T [H(y + kb - e) + H(y - kb + e)]. 

k= — oo 

Let us expand the function <&o(y) in a series of e. Applying the two-sided 
Laplace transform (?), one obtains: 

oo 

$(p,e)= J e-^y^(y,e) dy 

— OO 

Expanding the function 4>(p, e) in a series of e and performing the inverse 
Laplace transform, we obtain 

*o(y) = 2£$(y) + 2s Yl e"* (n) (y), (79) 

fc=l,3,5,... 
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CO 

where $>(y) = ^ S(y — kb). 

h— — oo 

Now let us consider the solution of the Eq. ([75)) in the form 

u = uq + eu\ + £ 2 ^2 + (80) 

Substituting ansatzes ([79]) , ([80)) into Eq. ([75]) and splitting the resulting 
equation with respect to e, we arrive at the recursive sequence of BVPs 

[l + 2£<l>(y)]u 0xx +u 0yy = 0, (81) 

<My), (82) 

Thus, in the zero approximation, we obtain the problem with one-dimensional 
fibers ([5T]h the influence of the width of the filaments is taken into account 
in the following approximations ([52]) . 

3 Summation of Asymptotic Series 

3.1 Analysis of Power Series 

We suppose that one obtains the following series as a result of an asymp- 
totic study: 

oo 

/(£)~^C n e n for e^O. (83) 

As it is known, the radius of convergence So of series ([53]) is determined by 
the distance to the nearest singularity of the function f(e) on the complex 
plane and can be found using the Cauchy-Hadamard formula: 

-=TST|C n | 1/n . 

If the nearest singularity lies on the positive real axis, then the coefficients 
C n have the usually one and the same algebraic sign, for example, 

— !— ~l+£ + £ 2 +£ 3 + .... 

1 - e 

If the nearest singularity is located on the negative axis, the algebraic signs 
of the coefficients C n are usually alternated, for example, 

- l- £ + £ 2 -£ 3 + .... 

1 + e 
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The pattern of signs is usually set pretty quickly. If there are several features 
of the same radius, that could happen to a real function with complex 
singularities necessarily occurring in complex conjugate pairs, then the rule 
of alternation of signs may be more complex, such as 



1 + 6 

Here we have a pattern of signs + H . To define eo may it is useful to 

apply the so-called Domb-Sykes plot (????). Let the function / have one 
of the nearest singularities at a point s = ±£q with an index of a, i.e. 



then we get 



(e ±e) a for 0,1,2,..., 

(e ± e) a ln(e ± e) for a = 0,1,2,..., 



C n 1 A 1 + a 



C n _i £o 



Constructing a graph of C n /C n -\ an the vertical axis and 1/n on the 
horizontal axis, one obtains the radius of convergence (as the reciprocal 
of the intercepts on the axis C n /C n _i), and then, knowing the slope, the 
required singularity. Figure 1 shows the numerical results for the function 

f(e) = e(l + e)(l + 2e)- 1 / 2 

9 3 „ 3 A 27 , 51 fi 191 7 359 o 

f — f -I <r<* p 4 -I p b -I <r' p» I 

+ 2 2 + 8 8 + 16 16 ' 

(84) 

starting with n = 7 points arranged in a linear relationship. 

If £o o r <^ are known from physical considerations, they can be used 
for the construction of the Domb-Sykes plot. If several singularities have 
the same convergence radius, so that the signs of the coefficients oscillate, 

1/2 

you can try to construct a dependence on the value (C n /C n _i) . If the 
radius of convergence tends to infinity and C n /C n -\ ~ k/n, then the ana- 
lyzed function has a factor exp(fce:), when C n /C n -i ~ k/n 1 ^ it has a factor 
exp(e p ). If the radius of convergence tends to zero, then the analyzed func- 
tion has an essential singularity and asymptotic expansion diverge. If the 
coefficients behave like C n -i/C n ~ l/(kn) then we can write C n ~ Ck n n\, 
where C is a constant. 

Knowledge of the singular solutions can eliminate them from the pertur- 
bation series and thus significantly improve its convergence. We describe 
some techniques for removing singularities. If the singularity lies on the 



18 



















. © - 














O 





























0.1 0.2 03 4 OS 1/n 



Figure 1. The Domb-Sykes plot for f(e) = e(l + + 2s;)- 1 / 2 . 



positive real axis, then it often means that the function f(s) is multivalued 
and that there is a maximum attainable point e = so- Then the inverse of 
the original function e = e(f) can be single valued. For example, consider 
the function 

f(e) = arcsine = e + ^ 3 + + -^e 7 + . . . , (85) 

the inverse of this function is 

e ~ / - -/ 3 + -20 f — f + . . . . (86) 

J 6 J 1 J 5040 ^ V J 

Numerical results are shown in Figure 2, where the solid line denotes the 
function arcsine, the dotted and dashed line shows the n-term expansions 
(|55j) and k- terms expansions ([86]) for different numbers of terms. It is evident 
that the expansion ([86]) allows a good description of the second branch of 
the original function. 
If 

/ ~ A(e - e) a for e -> e , < a < 1, 

the transition to the function J 1 / 2 removes the singularity. 
Consider the following example: 

e/2 FT" , 1 7 2 41 o 367 4 4849 5 

/(*) = e-/»>/l + 2i ~ 1 + -e - -/ + - £ 3 - — e 4 + — £ « + . . . . (87) 
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Figure 2. Application of the inversion of a power series. 



The radius of convergence of this expansion is equal to 1/2, while the 
radius of convergence of functions 



is infinite. 

Numerical results are shown in Figure 3, where the solid line denotes 
the function f(e) — e~ e l 2 \/\ + 2e, the dotted and dashed lines shows the 
n-term expansions (|57|) and the square roots of fc-term expansions ([88]) . 

In addition, knowing the singularity, one can construct a new function 
/m(^) (multiplicative extraction rule) 



The functions /m(^) and /a(£) should not contain singularities at eo. 
In many cases, one can effectively use the conformal transformation of the 
series, a fairly complete catalog of which is given in (?). In particular, it 
sometimes turns out to be a successful Euler transformation (???), based 
on the introduction of a new variable 




(88) 



f(s) = (e -s) a f M (s), 



or /a(£) (additive extraction rule) 



f(e) = A(e - e) a f A (e). 



e 



(89) 



£ 



1 - e/eo ' 
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0.5" I . . - , J i — , LI — . — . — . — I — 

o.o as to is £ 

Figure 3. An example of taking the root. 



Recast the function / in terms of e , / ~ ^d n e n , has the singularity 
pushed out at the point e = oo. For example, the function ([84]) is singular 
at the s = —1/2, which can be eliminated with the Euler transformation 
e = e/(l + 2s). The expansion of the function (|8^|) in terms of e is 

, !~ !~2 31 o 895 , 22591 , 

f (e) - 1 + -e + -e 2 e 3 e 4 e 5 + . . . . (90) 

M ; 2 8 48 384 3840 V y 

Some numerical results are shown in Figure 4, where the dotted and dashed 
line shows the n-term expansions and /c-terms in the expansion ([90]) . 
A natural generalization of Euler transformation is 

e 

£= (l-e/eo)"' 

where a is the real number. 

3.2 Pade Approximants and Continued Fractions 

; T/ie coefficients of the Taylor series in the aggregate have a lot more 
information about the values of features than its partial sums. It is only 
necessary to be able to retrieve it, and some of the ways to do this is to 
construct a Pade approximant" (?). Pade approximants (PA) allow us to 
implement among the most salient natural transformation of power series 
in a fractional rational function. Let us define a PA following ?. Suppose 
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we have the power series 

CO 

= (91) 

i=l 

Its PA can be written as the expression 

ao + aieH h a n e n 

f ^ {£) = l + 6ie+ ... + 6me m' ( 92 ) 

whose coefficients are determined from the condition 



• •+6 m e m )(co+ci£+c 2 £ 2 + . . . ) = a +a l£ +- ■ ■+a n e n +0(e m+n+1 ). 

(93) 

Equating the coefficients of the same powers e, one obtains a system of 
linear algebraic equations 

frraCn-ra+l + ^m-lC n -m+2 + C n +1 = 0; 
Dm — 1 — rn + 3 

+ c„ +2 = 0; 

: : : : (94) 

where Cj = for j < 0. 
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The coefficients ai can now be obtained from the Eqs. ([93]) by comparing 
the coefficients of the powers e: 

a 
ai 



where p = min(n, m). Eqs. ([93|) , ([94]) are called Pade equations. In the case 
where the system (|9l|) is solvable, one can obtain the Pade coefficients of 
the numerator and denominator of the PA. Functions f[ n / m ] ( e ) a ^ different 
values of n and m form a set, which is usually written in the form of a table, 
called the Pade table (Table 3): 



m 





1 


2 







f\o/o] (e) 


f\i/o] (e) 


f\2/0](. £ ) 




1 


ijo/nO) 


/fi/n(^) 


f\2/l](e) 




2 


f\o/2] (e) 


f\i/2] (e) 

















Table 3. Pade table. 

The terms of the first row of the Pade table correspond to the finite sums 
of the Maclaurin series. In case of n = m one obtains the diagonal PA, the 
most common practice. Note that the Pade table can have gaps for those 
indices n, m, for which the PA does not exist. We note some properties of 
the PA (???): 

1. If the PA at the chosen m and n exists, then it is unique. 

2. If the PA sequence converges to some function, the roots of its denom- 
inator tend to the poles of the function. This allows for a sufficiently 
large number of terms to determine the pole, and then perform an 
analytical continuation. 

3. The PA has meromorphic continuation of a given power series func- 
tions. 

4. The PA on the inverse function is treated the PA function inverse 
itself. This property is called duality and more exactly formulated as 
follows. Let 

q(s) = f-\s) and /(0) + 0, then q [n/m] (s) = f^ /m] (e) , (96) 



2 = 1 



(95) 
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provided that one of these approximations there. 

5. Diagonal PA are invariant under fractional linear transformations of 
the argument. Suppose that the function is given by their expansion 
([9T]). Consider the linear fractional transformation that preserves the 
origin W = jt^, an d the function q(W) = f(s). Then <?[ n / n ], pro- 
vided that one of these approximations exist. In particular, the diag- 
onal PA is invariant concerning Euler transformation ([89]) . 

6. Diagonal PA are invariant under fractional linear transformations of 
functions. Let us analyse a function (|9T|) . Let 



If c + df{e) ^ 0, then 

Q[n/n] 0=0 



bf[n/n] (g) 

' df[ n / n ] (e) ' 



provided that there is f[ n / n ] (e). Because of this property infinite values 
of PA can be considered on a par with the end. 

7. The PA can get the upper and lower bounds for f[ n / n ](e). For the 
diagonal PA one has the estimate 

/[n/n-l](£) < f[n/n}(e) < /[n/n+1] (^) • (97) 

Typically, this estimate is valid for the function itself, i.e. f[ n / n ] (e) 
in Eq. ([97|) can be replaced by f(e) . 

8. Diagonal and close to them a sequence of PA often possess the prop- 
erty of autocorrection (??). It consists of the following. To determine 
the coefficients of the numerator and denominator of PA have to solve 
systems of linear algebraic equations. This is an ill-posed procedure, 
so the coefficients of PA can be determined with large errors. How- 
ever, these errors are in a certain sense of self-consistency, the PA can 
approximate the searching function with a higher accucary. This is a 
radical difference between the PA and the Taylor series. 

Autocorrection property is verified for a number of special functions. At 
the same time, even for elliptic functions the so-called Froissart doublets 
phenomenon arises, consisting of closely spaced zeros and poles to each 
other (but different and obviously irreducible) in the PA. This phenomenon 
is not of numerical nature, but due to the nature of the elliptic function (?). 
Thus, in general, having no information about the location of the poles of 
the PA, but relying solely on the PA (computed exactly as you wish), we 
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can not say that you have found a good approximation for the approximated 
function. 

To overcome these defects several methods are suggested, in particular, 
the smoothing method (?). Its essence is that instead of the usual-term 
diagonal PA for complex functions f[ n / n ](e) = Pn(^)/Qn(^) the following 
expression is used 

- / x _ q n (s)p n (s) + q n -i(s)p n -i(s) 

J[n/n]\ e ) — -r-- ( v -j-r- , v 

Here / denotes complex conjugation of /. 

Now consider the question: in what sense can the available mathemati- 
cal results on the convergence of the PA facilitate the solution of practical 
problems? Gonchar's theorem (Gonchar, 1986) states: if none of the di- 
agonal PA f[ n /n]{e) has poles in the circle of radius R, then the sequence 
f[n/n]( £ ) is uniformly convergent in the circle to the original function f(s). 
Moreover, the absence of poles of the sequence of the f[ n / n ](e) in a circle of 
radius R must be original and confirm convergence of the Taylor series in 
the circle. Since for the diagonal PA invariant under fractional linear maps 
we have s — » ^q^, the theorem is true for any open circle containing the 
point of decomposition, and for any area, which is the union of these circles. 
A significant drawback in practice is the need to check all diagonal PA. For 
example, the following theorem holds (?): Suppose the sequence of diago- 
nal Pade approximants of the function w(e), which is holomorphic in the 
unit disc and has no poles outside this circle. Then this sequence converges 
uniformly to w(e) in the disc \z\ < ro, where 0.583R < ro < 0.584R. 

How can we use these results? Suppose that there are a few terms of the 
perturbation series and someone wants to estimate its radius of convergence 
R. Consider the interval [0, £o], where the truncated perturbation series and 
the diagonal PA of the maximal possible order differ by no more than 5% 
(adopted in the engineering accuracy of the calculations). If none of the 
previous diagonal PA does not have poles in a circle of radius £q> then it is 
a high level of confidence to assert that R>£q- 

The procedure of constructing the PA is much less labor intensive than 
the construction of higher approximations of the perturbation theory. The 
PA is not limited to power series, but to the series of orthogonal polynomi- 
als. PA is locally the best rational approximation of a given power series. 
They are constructed directly on its odds and allow the efficient analytic 
continuation of the series outside its circle of convergence, and their poles 
in a certain sense localize the singular points (including the poles and their 
multiplicities) of the continuation function at the corresponding region of 
convergence and on its boundary. This PA fundamentally different from 
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rational approximations to (fully or partially) fixed poles, including those 
from the polynomial approximation, in which case all the poles are fixed 
in one, infinity, the point. Currently, the PA method is one of the most 
promising non-linear methods of summation of power series and the local- 
ization of its singular points. Including the reason why the theory of the 
PA turned into a completely independent section of approximation theory, 
and these approximations have found a variety of applications both directly 
in the theory of rational approximations, and in perturbation theory. Thus, 
the main advantages of PA compared with the Taylor series as follows: 

1. Typically, the rate of convergence of rational approximations greatly 
exceeds the rate of convergence of polynomial approximation. For ex- 
ample, the function e £ in the circle of convergence approximated by 
rational polynomials P n (s)/Q n (e) in 4 n times better than an algebraic 
polynomial of degree 2n. More tangible it is property for functions 
of limited smoothness. Thus, the function \e\ on the interval [—1,1] 
can not be approximated by algebraic polynomials, so that the or- 
der of approximation was better than 1/n, where n is the degree of 
polynomial. PA gives the rate of convergence ~ exp (— y/2nj. 

2. Typically, the radius of convergence of rational approximation is a 
large compared with power series. Thus, for the function arctan(x) 
Taylor polynomials converge only if \e\ < 1, and AP - everywhere in 
C \ ((—zoo, —i] U [i, ioc)). 

3. PA can establish the position of singularities of the function. 
Similarly, the PA method is a method of continued fractions (?). There are 
several types of continued fractions. The regular C-fraction has the form of 
an infinite sequence, in which N-th term can be written as follows 

h{e) = a + 2Lj (98) 

l + 

1 + c N e 

The coefficients q are obtained after the decomposition of expression 
([98]) in a Maclaurin series and then equating the coefficients of equal powers 
of e. When a = one obtains the fraction of Stieltjes or S- fraction. For the 
function of Stieltjes 

o 
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the coefficients of expansion ([98]) have the form: a = 0, cq = 1, C2 n -i = 
c 2n = n, n > 1. 

Description of the so-called J-, T-, P-, R-, g- fractions, algorithms for 
their construction and the range of applicability are described in detail in 

Continued fractions are a special case of continuous functional approxi- 
mation (?). This is the sequence in which the (n + l)-th term c n (e) has the 
form n-th iteration of a function F(e). For the Taylor series one obtains 
F(e) = 1 + £, for the continuous fraction F(e) = j^. If F(e) = exp(£) one 
obtains a continuous exponential approximation 

c n (e) = clq exp {ai£:exp [a2£ • • . exp (a n e)}} , 



for F(e) = VTTe 



c n (e) = a Y 1 + ai£\l 1 + a 2 e\/ 1 + . . . a n -is^l + a n e, 
for F(e) = In 1 + e 

c n (e) = do In {a\£ In \a 2 e ... In (a n e)]} . 

In some cases, such approximations can converge significantly faster than 
power series. 

As an example, we note the solution of the transcendental equation 

x = e In x 

for large values (?, §3.4.9)) A: 

xo = elne; x\ = e ln(elne); x 2 = e ln[e ln(e lne)]; 

4 Some Applications of Pade Approximants 

4.1 Accelerating Convergence of Iterative Processes 

The efficiency of PA or other methods of summation depends largely on 
the availability of higher approximations of the asymptotic process. Some- 
times they can be obtained by using computer algorithms (?), but in general 
it remains an open question. Iterative methods are significantly easier to 
implement. As a result of an iterative procedure a sequence of S n is ob- 
tained. Suppose that it converges and has the limit value. We introduce 
the parameter a obtained by the ratio 

i • $n+l — $>n 

a= lim — o q • 

n— )-oo Or, — O 
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It's called superlinear convergence, if a = 0, a linear for a < 1 and 
logarithmic at a = 1. The biggest issues are, of course, logarithmically 
convergent sequences - alas, widespread practice. Very often linearly con- 
vergent sequences are also a problem. Therefore, it is often necessary to 
improve the convergence. One method of improving the convergence is to 
move to a new sequence T n with the aid of a transformation so that 



run 



>oo S n — S 



In such cases we say that the sequence T n converges faster than sequence 
S n . There are linear and nonlinear methods to improve convergence. Linear 
methods are described by formulas 



T n — y a n iSi, n — 0, 1, 2, . . . 



where the coefficients a n { do not depend on the terms of the sequence S n 
constant. 

Since linear methods improve the convergence of a restricted class of 
sequences, currently nonlinear methods the most popular ones. Among 
them thr Aitken method (?) stands out for its easiness, which described by 
the formula 



T n = S n - {Sn ^ S " )(S " n = 0,l,2,.... (99) 



The Aitken method accelerates the convergence of all linear and many 
of logarithmically convergent sequences. It is very easy to calculate, and 
in some cases it can be applied iteratively. A natural generalization of the 
Aitken transformation is the Shanks transformation (?) 



T; h = -gy, (100) 
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where 



D 



(i) 

kp 



D 



(2) 
kp 



Sp—k 










k 


A^p-fe+i . 


.. AS P 


AS P - 


i 


AS P 


• • A^+zc-i 


1 




1 


1 


AV 


k 


AaS'p-zc+i 


AS P 


AS P - 


i 


AS P 


• • AS p +k-i 



ASk = /Sfc+i - Sk 

Eq. (|100p is called the Shanks transformation of the order k of the 
sequences to the sequence Sk- For k = 1 one obtains the Aitken transform 
([99]) . Shanks method requires the calculation of determinants, which is not 
always easy. One can use also Wynn algorithm, which described by the 
formulas: 



0, Tr 



(n) 



Sni Tn 



(n) 



fc+1 



l k+l 



1 



i(n+l) 



(n)- 



(101) 



The Wynn algorithm is related to the transformation of Shanks ()100j) in the 
following way: 



r (n) 



T fc W (S n ), Tr 



(n) 
2/c+l 



1 



T, 



(A5 n ). 



The Wynn algorithm is a quadratic convergent method for solving sys- 
tems of nonlinear equations (?). There are many other techniques for ac- 
celerating of sequences' convergence. One can use them consistently, for 
example, to convert the original sequence into a linearly convergent one, 
and then apply the method of Aitken. One can also use different methods 
to improve convergence, at each stage by comparing the obtained results 
(?). All the described methods have a close relationship with the PA. The 
Aitken method corresponds to the PA [n/1] , the Shanks method to the PA 
\p/k], for the method of Wynn one obtains = [n + k/k). 



4.2 Removing Singularities and Reducing the Gibbs Effect 

Consider the problem of uniform plane flow of an incompressible inviscid 
fluid pasting a thin elliptic airfoil (\x\ < 1, \y\ < e, e <C 1). The expression 
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for the relative velocity q* of the flow is such (?, §4.4): 

= ^)^)^\ (102 ) 

where V is the free-stream speed. 

The splitting of the r.h.s. of Eq. (|102[) in a series of £ can be expressed 

as: 

9*(* >e ) = 1 + e - \e 2 ^ 2 ~ \^Y^ + ■ ■ ■ ■ ( 103) 

This expression diverges at x = 1. No wonder it is not: the expansion 
(|103p is obtained as a result of the limiting process lim q(x,e), x > 1, and 

to get the value of it is necessary to perform the limit as lim 

for e > 0. Divergence of series (|103p when x = 1 indicates that the limit 
processes cannot be interchanged. Now let us apply PA to the r.h.s. of Eq. 
([103)) and then pass to the limit x —> 1. After trying various options, we 
conclude that the best result is given by the PA 

= <lz£0±£>. (10 4) 

Numerical results for e = 0.5 are shown in Figure 5, where the dashed 
line denotes the solution (|103p , curves 1 and 2 the exact solution (|104j) 
and the PA (|104p . It is seen that the use of PA significantly improves the 
accuracy of the approximate solution. 

PA can be also successfully applied for the suppression of the Gibbs 
phenomenon (????) . Consider, for example, the function sign x: 



sign x 



-1, -7r < x < 0, 

1, < X < 7T. 

Its Fourier series expansion has the form 

4^sin(2j + l)x 
sign x = - > — y — J — (105) 

& 7T^ 27 + 1 V ) 

3=0 J 

Direct summation of series (|105p leads to the Gibbs effect in the neigh- 
borhood of x — 0, while the defect of convergence reaches 18%, i.e. instead 
of 1 one obtains the value of 1.1789797. . .. Diagonal PA for series (|105j) 
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Figure 5. Removing singularities by the PA. 



can be written as follows: 



sign x [N/N] 



N-l/2 

g 2 j + isin((2j + l)x) 

3=0 

N/2 

1 + iZ s 2j cos(2jx) 

j=0 



where 



q 2j+1 = -(2j + l) 

7T 



[N/2] 
E 



S2i 



S2t 



(106) 



(2j + l)2 ,tl (2j + 1)2 - (2i) 2 

(AH) 4 (2A/" + 2z)!(2A"-2i)! 
(AT - 1)!(JV + 1)!(JV - 2z)!(AT + 2z)![(2JV)!] 2 ' 



Numerical studies show that the Gibbs effect for PA (|106p does not 
exceed 2% (?). 



4.3 Localized Solutions 

We consider the stationary Schrodinger equation 

V 2 0, y) - u(x, y) + ?i 3 0, 2/) = 0. 



(107) 
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We seek the real, localized axisymmetric solutions of the Eq. (|107j) . In polar 
coordinates (£, 0) we construct a solution ip(£) which does not depend on 0. 
As a result, we obtain the BVP 



<p"(0 + V(0-v(0 + v 3 (0 



(108) 



)>co 



lim <£>(£) 



0, 
0. 



(109) 
(110) 



BVP (|108|) - (|110p can be regarded as an eigenvalue problem, and the role 
of an eigenvalue is an unknown quantity A = (p(0). This problem plays 
an important role in nonlinear optics, quantum field theory, the theory of 
magnetic media. As shown in (?, p. 12-16), BVP (|108|) - (|110|) has a countable 
set of "eigenvalues" A n , the solution A n ) has exactly n zeros, and the 
solution (p(£,Ao) has no zeros and decreases monotonically on £. That is 
the last solution, which is most interesting from the standpoint of physical 
applications, and we will concentrate on obtaining it. 

The problem of computing the decaying solutions of BVP (jl08[) - ([110p 
is identical to the problem of computing homoclinic orbits in the 3D phase 
space for the nonlinear oscillator, or equivalent ly, for computing the initial 
conditions for these orbits. 

Since the sought solutions are expected to be analytical functions of £, 
they can be expressed in Maclaurin series about £ = 0: 



Substituting ansatz (|111|) into Eq. (|108|) , producing a splitting of the 
powers of the £ and solving the relevant equations, one obtains (?): 



C 2 = 0.25Ao(l - Al); C7 4 = 0.25<7 2 (7; C 6 = -± £— ; C 8 = 



CO 




(111) 



^(CC 6 - QA C 2 C 4 - C|); 
Cio = -O.Ol^Cg + QA C 2 C 6 



3AqC1 + 3C2C4); 



C12 = -^(-CCio + 6A C 2 C 8 + 6A C 4 C 6 + 3(7 2 ci), 
where C = 1-3AZ. 
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Then we construct PA for the truncated series (jllll) 

JV 

v(0 = £5 (ii2) 

i + e b 2j e k 

k=i 

All coefficients in Eq. (|112[) can be parameterized in terms of Aq, a 2 j = 
a2j(Ao), = b2j(Ao) and the PA becomes a one-parameter family of 
analytical approximations of the searching solution. Then we compute the 
value of Aq for which PA (|112p decays to zero as £ tends to infinity. It gives 
us condition 

a 2i (A ) =0, b 2j (A ) ^ 0, (113) 

One can compute the PA (|112|) , then imposed the condition (|113[) and 
obtained the following convergent values of Aq for varying orders 2 AT: 



TV 


1 


2 


3 


4 


A 


±v / 3 


±2.20701 


±2.21121 


±2.21200 



The numerical solution gives Aq « ±2.206208, the difference between 
numerical and analytical solutions for N = 4 is 0.26%. 

4.4 Hermit e-Pade Approximations and Bifurcation Problem 

PA can successfully work with functions having poles. However, it often 
becomes necessary to explore functions with branch points, and construct 
all their branches. In that case, one can use Hermite-Pade approximations 
(?). Suppose it comes to a function with the expansion 

oo 

/(£) = $>„£", (114) 
n=l 

and we managed to find the first few coefficients of this series 

N 

fN(e) = ^u n £ n . 

n=l 

If it is known that this function has a branch point, we can try to transform 
the original series (|1 14j) in an implicit function 

F(e,f)=0, 

and determine all required branches of it. 
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For this purpose we construct a polynomial F p (e, f) of degree p > 2 

p m 
m=l k=0 

It was assumed Co.i = 1, and the remaining coefficients must be deter- 
mined from the condition 

F p (e, /*(£)) = 0(s N+1 ) at £ -> 0. (115) 

Polynomial F p contains 0.5 (p 2 + 3p — 2) unknowns, the condition (|115|) 
yields TV linear algebraic equations, therefore, should A/" = 0.5 (p 2 + 3p — 2). 
Once the polynomial F p is found, one can easily find p branches of the 
solution from the equation 

F p = 0. 

For the analysis of bifurcations of these solutions one can use Newton's 
polygon (?). If a priori information about the searching function is known, 
it can be taken into account for constructing the polynomial F p . 

4.5 Estimates of Effective Characteristics of Composite Materi- 
als 

We consider a macroscopically isotropic 2D composite material consist- 
ing of a matrix with inclusions. The aim is to determine the effective con- 
ductivity q from the known matrix (qi) and inclusions (q 2 ) conductivities 
and the volume fraction ((f). As shown in (?), if we take e = ^ — 1 as a 
small parameter value, the required effective conductivity can be written as 
follows 

— = l + y>e- 0.5(^(1 - (p)e 2 + 0(s 3 ). (116) 
Qi 

Using the first two terms of expansion (|116|L one obtains 

t~~ < - < 1 + <pe, 
l + pe qi 

hence the bounds of Wiener 



+ — ) <q<(l- (p)qi + <pq 2 . 
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3 2 n I n n+] 

Figure 6. A chain of elastically coupled masses. 



4.6 Continualization 

We study a chain of n + 2 material points with the same masses m, 
located in equilibrium states in the points of the axis x with coordinates 
jh(j = 0, 1, . . . , n, n + 1) and suspended by elastic couplings of stiffness c 
(Figure 6) (?). 

Owing to the Hooke's law the elastic force acting on the j-th mass is as 
follows 



= c[Vj+i(t) -Vj(t)] -c[yj(t) -yj-i{t)] = c[y j - 1 (t)-2y j (t)+y j+1 (t)], 

where j = 1, 2, . . . , n and yj(t) is the displacement of the j-th material point 
from its static equilibrium position. 

Applying the second Newton's law one obtains the following system of 
ODEs governing chain dynamics 

ma jtt (t) = c((Tj+i -2(jj +cr i _i), j = 1,2, ...,n. (117) 

Let us suppose the following BCs 

<ro(t) = * n+1 (t) = 0. (118) 

For large values of n usually continuum approximation of discrete problem 
is applied. In our case it takes the form of 

ma t t(x,i) = ch 2 a xx (x,t), (119) 
(j(0, t) = a(£,t) = 0. (120) 

Formally, one can rewrite Eq. (|117|) as pseudo-differential equation: 

m w +4csm -y^r = - (121) 
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Pseudo-differential operator can be split into the Maclaurin series as follows 

. o / ih d \ 1 „ h 2k d 2k 

sm — T«Z = -o E 



2 dx) 2 fc f i (2fc)! &r 2fc 



h 2 d 2 ( h 2 d 2 /i 4 a 4 /i 6 a 6 
i + — ■ 



4 &c 2 V 12 dx 2 360 5x 4 10080 5a; 6 

(122) 

With only keeping the first term in the last line of Eq. (|122p , one obtains 
a continuous approximation ([119[) . Keeping the first three terms in Eq. 
(|122p , the following model is obtained 

d 2 a _ u2 ( d 2 h 2 d 4 h 4 d 6 \ 
m dt 2 ~ \dx 2 + 12 dx 4 + 360 dx* ) ( 3) 

In the case of periodic BCs for a discrete chain one obtains the following 
BCs for Eq. ([125)1 : 

cr = (J xx = (?xxxx = for x = 0, £. (124) 

BVP (|123p . (|124|) is of the 6th order in spatial variable. Using PA we can 
obtain a modified continuous approximation of the 2nd order. If only two 
terms are left in in the last line of Eq. (|122p , then the PA can be cast into 
the following form: 

d 2 

d 2 h 2 d 4 ^ 
dx 2 12 dx 4 ~ h 2 d 2 ' 
1 ~12dx~ 2 

For justification of this procedure Fourier or Laplace transforms can be 
used. 

The corresponding so-called quasicontinuum model reads 



12 dx 2 ) 



(j u - ch 2 a xx = (125) 
The BCs for Eq. (p5l) have the form (IT2Q1) . 



4.7 Rational interpolation 

Here we follow (?). A simple way to approximate a function is to choose 
a sequence of points 

a = xq < x\ < X2 • • • < x n = 6, 
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and to construct the interpolating polynomial p n (x) 

Pn(xi) = f(xi), i = 0,1,2,..., n. 

However, as is well-known p n (x) may not be a good approximation to /, 
and for large n ^> 1 it can exhibit wild oscillations. If we are free to choose 
the distribution of the interpolation points one remedy is to cluster them 
near the end-points of the interval [a, 6], for example using various kinds of 
Chebyshev points. 

Sowe have to make do with them, and then we need to look for other 
kinds of interpolants. A very popular alternative nowadays is to use splines 
(piecewise polynomials), which have become a standard tool for many kinds 
of interpolation and approximation algorithms, and for geometric modeling. 
However, it has been known for a long time that the use of rational functions 
can also lead to much better approximations than ordinary polynomials. In 
fact, both, polynomial and rational interpolation, can exhibit exponential 
convergence when approximating analytic functions. In "classical" rational 
interpolation, one chooses some M and N such that M + N = n and fits 
a rational function of the form to the values where pm and qw 

are polynomials of degrees M and N respectively. If n is even, it is typical 
to set M + N = 77, and some authors have reported excellent results. The 
main drawback, though, is that there is no control over the occurrence of 
poles in the interval of interpolation. 

? suggested that it might be possible to avoid poles by using rational 
functions of higher degree. They considered algorithms which fit rational 
functions whose numerator and denominator degrees can both be as high 
as n. This is a convenient class of rational interpolants because every such 
interpolant can be written in so-called barycentric form 



r(x) 



n \. 

E -^-f(xi) 

i=0 X Xi 

E -z-fixi) 

i—0 X Xi 



for some real values A^. Thus it suffices to choose the weights Xi in order 
to specify r, and the idea is to search for weights which give interpolants r 
that have no poles and preferably good approximation properties. Various 
approach is described in ?, in particular, one can choose A^ = (— 1)*, i = 
0,1,2,. 

4.8 Some Other Applications 

PA is widely used for the construction of solitons and other localized 
solutions of nonlinear problems, even in connection with the appeared term 
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"padeon" (??). As a simple model, we consider the nonlinear BVP 

y"-y + 2y 3 = 0, (126) 
y(0) = 1, y(oo)=0, (127) 

that has an exact localized solution 

y = cosh~ 1 (x) (128) 

Quasilinear asymptotics give a solution in the following form: 

y = Ce~ x (1 - 0.25C 2 e~ 2x + 0.0625C 4 e - 4 * + ...), C = const. (129) 



It is easy to verify that with reconstructing the truncated series ([129)) in the 
PA and to determine the constant C from the BCs (|127|) and we arrive at 
the exact solution (|128[) . 

It is also interesting to use the PA to the problems with the phenomenon 
of "blow-up", when the solution goes to infinity at a finite value of the 
argument. For example, the Cauchy problem 

^-=ax + ex 2 , x(0) = l, 0<e«a«l, (130) 
has the exact solution 

X(t) = aeMat \ r , (131) 

a + e — eex.-p(at) 

which tends to infinity for t — >• ln[(a + s)/e\. 
Regular asymptotic expansion 

x(t) ~ exp(at) — soT 1 exp(at)[l — exp(at)] + . . . 

can not describe this phenomenon, but the use of the PA gives the exact 
solution (dSH) . 

PA allows to expand the scope of the known approximate methods. For 
example, in the method of harmonic balance the representation of the solu- 
tion of a rational function of the type 

N 

{A n cos[(2n + l)ut] + B n sin[(2n + 

= 

1+5^ {Cm cos(2mut) + D m sin(2mcjt)} 

n=0 
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substantially increases the accuracy of approximation (??). PA can be used 
effectively to solve ill-posed problems. This could include reconstruction 
of functions in the presence of noise (??), various problems of dehomoge- 
nization (i.e., determining the components of a composite material on its 
homogenized characteristics) (?), etc. We must also mention 2D PA (?). 



5 Matching of Limiting Asymptotic Expansions 

5.1 Method of Asymptotically Equivalent Functions for Inver- 
sion of Laplace Transform 

This method was originally proposed by Slepian and Yakovlev for the 
treatment of integral transformations. Here is a description of this method, 
following ?. Suppose that the Laplace transform of a function of a real 
variable f(t) is: 



F(s) = jf(t)e 





To obtain an approximate expression for the inverse transform, it is 
necessary to clarify the behavior of the transform the vicinity of the points 
5 = and s = oo and determine the nature and location of its singular 
points are on the exact boundary of the regularity or near it. Then the 
transform F(s) replaced by the function Fq(s), allowing the exact inversion 
and satisfying the following conditions: 

1. Functions Fq(s) and F(s) are asymptotically equivalent at s — » oo and 
8—^0, i.e. 

Fq(s) ~ F(s) at s — ^ and s — ^ oo. 

2. Singular points of the functions Fq(s) and F(s), located on the exact 
boundary of the regularity, coincide. 

The free parameters of the function Fq(s) are chosen in such a way that 
they satisfy the conditions of the approximate approximation of F(s) in the 
sense of minimum relative error for all real values s > 0: 



mm < max 



F (s,ai,a 2 , • • • _ 1 



F(s) 



(132) 



Condition (|132p is achieved by variation of free the parameters ol{. Often 
the implementation of equalities 



oo oo 

J F (s) ds = J F(s) ds 





39 



or Fq ~ Fq at s —> leads to a rather precise fulfillment of the requirements 
(|132p . Constructed in such a way the function is called asymptotically 
equivalent function (AEF). 

Here is an example of constructing AEF. Find the inverse transform if 
the Laplace transform is the modified Bessel function (?, Sect. 9): 



K (s) = -ln(s/2)I (s) 



00 2k 

V — 

2 2fc (fe!) 2 



tf(Jfe + l) 



(133) 



where *&(z) is the psi (digamma) function (?, Sect. 6). 

For pure imaginary values of the argument s(s = iy\ < \y\ < 00) 
function Kq(s) has no singular points. Consequently, we can restrict the 



study of its behavior for s —> and s 
expressions are (?, Sect. 9): 



00. The corresponding asymptotic 



K (s) 
Ms) 



In — h 7 
. 2 . 



■0(a), 



(134) 



OO, 



where 7 is the Euler's constant (7 = 1, 781 . . . ) (note the typo in the first 
formula (1131) in ?). 

The analyzed Laplace transform has a branch point of the logarithmic 
type, branch point of an algebraic type, and an essential singularity. These 
singular points need to be stored in the structure of the zero approximation. 
The most simple way to obtain such a structure, combining two asymptotic 
representations (|134p so that they are mutually do not distort each other 
and contain free parameters, which could be disposed of in the future. As 
a result, we arrive at the zero approximation 



Fo(s) 



In- 



a 



1 



2 v^T^J 



(135) 



where a and j3 are the free parameters. 

It is easy to see that expression (|135p has the correct asymptotic be- 
havior s — » 00. The free parameters are determined from the condition of 
coincidence of the asymptotics of the functions Ko(s) and Fo(s) for s —> 
and the equality of integrals 



00 00 
J F (s) ds = J K (s) ds. 
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As a result of calculations one obtains a system of transcendental equa- 
tions 



irm+ vi = ln2_7; 

lna-e a Ei(-a)+ 7 +^e^ [l-erf(>/2)] = J, 
v 2 2 

where Ei(. . . ) is the the sine integral (?, Sect. 5), erf(. . . ) is the the error 
function (?, Sect. 7) (note typo in these formulas in ?). 

Solving the prescription system numerically, one finds a = 0.3192, f3 = 
0.9927. 

Then the approximate inverse transform can be written as follows: 

fM = { 1 - - 1)1 + }»(«-»)• U36, 

The exact expression for the function /(£) is: 

/(t) = ^=jH(t - 1) (137) 

Comparison of exact (|137p (solid line) and approximate (|136p (dotted line 
with circles) inversions is shown in Figure 7. As it can be seen, a satisfactory 
result is obtained even in the zero approximation. Analogously, one can 
construct AEFs for inverse sine and cosine Fourier transforms, Hankel and 
other integral transforms. 



5.2 Two-point Pade Approximants 

The analysis of numerous examples confirms: usually implemented a sort 
of "complementarity principle": if for e —> one can construct a physically 
meaningful asymptotics, there is a nontrivial asymptotics and e — »• oo. The 
most difficult in terms of the asymptotic approach is the intermediate case 
of e ~ 1. In this domain numerical methods typically work well, however, if 
the task is to investigate the solution depending on the parameter £, then it 
is inconvenient to use different solutions in different areas. Construction of 
a unified solution on the basis of limiting asymptotics is not a trivial task, 
which can be summarized as follows: we know the behavior of functions in 
zones I and III (Figure 8), we need to construct it in the zone II. For this 
purpose one can use a two-point Pade approximants (TPPA). We give the 
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f(t) 

5 T" 



10 



Figure 7. Comparison of the exact Laplace transform inversion with the 
treatment by the method of AEFs. 



definition following ?. Let 

oo 

F(e) = ^2 c i e% at £ ~^ °' 

i=0 

oo 

F(e) = CjS~ l at e 



oo. 



i=0 



(138) 
(139) 



Its TPPA is a rational function of the form 

f , v _ ap + aig j h a n g n 

/[n/m](£j " 1 + he + • • • + b m e™ 

with /c < m + n — 1 coefficients which are determined from the condition 

(1 + he H h frm£ m ) (co + cie + c 2 £ 2 + . . . ) = a + ai£ H h a n £ n 

and the remaining m + n — fc coefficients of a similar condition for e -1 . 

As an example of TPPA using for matching of limiting asymptotics, 
consider the solution of the van der Pol equation: 

x + ex (x 2 — l) + x = 0. 
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Asymptotic expressions of the oscillation period for small and large val- 
ues of e are (?): 

/ p2 t 4 \ 

T = 2tt 1 + — - — — at e^O, (140) 



16 3072 , 

T = e(3-21n2) at e oo. (141) 

For constructing TPPA we use the four conditions at £ — and the two 
conditions at e —> oo, then 

^/ \ ap + gig + a 2 g 2 + a3g 3 , 1/10 x 
T(£)= l + M + ' (M2) 

where 

7r 2 (3-21n2) 7r(3-21n2) 2 
ao = 27r; a\—— — — — —5 5-; a2 — 



4(3-21n2) 2 -7T 2 ' " 2[4(3-21n2) 2 -7r 2 ] 
7r 2 (3-21n2) 7r(3-21n2) 

a 3 = 77 r^/o 01 r,x9 onS "I = 



16[4(3-21n2) 2 -7r 2 ]' 1 2 [4(3 - 2 In 2) 2 - tt 2 ] ' 
2 16 [4(3 -2 In 2) 2 -tt 2 ]' 
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e 


T numerical 


T numerical 


1 


6.66 


6.61 


2 


7.63 


7.37 


3 


8.86 


8.40 


4 


10.20 


9.55 


5 


11.61 


10.81 


6 


13.06 


12.15 


7 


14.54 


13.54 


8 


16.04 


14.96 


9 


17.55 


16.42 


10 


19.08 


17.89 


20 


34.68 


33.30 


30 


50.54 


49.13 


40 


66.50 


65.10 


50 


82.51 


81.14 


60 


98.54 


97.20 


70 


114.60 


113.29 


80 


130.67 


129.40 


90 


146.75 


145.49 


100 


162.84 


161.61 



Table 4. Comparison of numerical results and calculations using the TPPA. 



Table 4 shows the results of the comparison of numerical values of the 
period, given in ?, with the results calculated by formula (|142|) . 

Now we construct inverse Laplace transform with the TPPA. Let the 
original function is as follows: 

/(*) = (l+t 2 )- - 5 . (143) 

Asymptotics of this function looks like: 

su\c±f l-0.5t 2 + ... at t->0, 
t~ 1 + ... at t^oo. 

TPPA in this case can be written as: 

.... 1 + 0.5* ..... 

f{t) = 1 + 0.5* + 0.5^ (144) 

Numerical results are shown in Figure 9. An approximate inversion (|144|) 
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f(t) 

0.8 
0.6 
0.4 
0.2 












































2 A 
Figure 9. Exact 


\ 6 t 8 

and approximate 



(upper curve) agrees well with the original (|143p (lower curve) for all values 
of the argument. 

Other example on the effective use of the TPPA see ????? . 

5.3 Other Methods of AEFs Constructiong 

Unfortunately, the situations where both asymptotic limits have the form 
of power expansions are rarely encountered in practice, so we have to resort 
to other methods of AEFs constructing. Consider, for example, the BVP 



ey xx -xy = ey, y(0) = 1, y(oo) = 0,£ «1. 

Solution for small values of x can be written as follows: 

1 



(145) 



(146) 



where £ = xe 3 ? a is the arbitrary constant. 

The solution for large values of x is constructed using the WKB method 

(?) 



y = b£ 4 exp 



3^ 



(147) 



where 6 is an arbitrary constant. 

Now we match these asymptotics. Because of the exponential in Eq. 
(|147p using TPPA in the original form is not possible. Therefore, we con- 
struct AEF, based on the following considerations: for large values of the 
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variable £ exponent from Eq. (|147|) is taken into account in its original 
form, and for small values of the variable £ it is expanded in a Maclaurin 
series. Constructed in this way AEF has the form 

2 .3 2.5 32 , 

y a = ^ 32. 3 5 exp (--^) (148) 

5 6 



The coefficients a and 6 in Eq. (|148j) still remain uncertain. For calculation 
of these constants one can use some integral relations, for example, obtained 
from Eq. (|148j) by multiplying them with the weighting functions 1, x, x 2 , . . . 
and further integration over the interval [0, oo). In the end, such values of 
the constants are found: 



a = ^#' b = ^h (149) 

V 

Numerical calculations show that the formula (|148p with constants (|149p 
approximates the desired solution in the whole interval [0, oo) with an error 
not exceeding 1.5%. 

When choosing the constants one can use other methods, then a lot 
depends on the skill of the researcher. Of course, it is necessary to ensure 
the correct qualitative behavior of AEFs, avoiding, for example, do not 
correspond to the problem of zeros of the denominator. To do this, one 
can vary the number of terms in the asymptotics and the numerator and 
denominator constructed uniformly suitable solutions. In general form the 
method of rational AEF can be described as follows (?). Let us assume that 
function f(z) has the following asymptotics: 

f(z) = F(z) at z -> oo, (150) 

and 

oo 

f(z) = J2c i z i at z^O. (151) 

i=0 

Then the AEF can be produced from the Eqs. (|152p , (|153p as follows: 

m 

E ot i (z)z i 

/(*)»^ at ( 152 ) 

i=0 
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where c^, ft are considered not as constants but as some functions of z. 
Functions c^(z) and P%{z) are chosen in such a way that: 

1. the expansion of AEF (|152p in powers of z for z — >• matches the 
perturbation expansion (|151|) : 

2. the asymptotic behavior of AEF (|152p for z — >> oo coincides with the 
function F(z) ([T50D . 

In the construction of AEFs a priori qualitative information is very im- 
portant. For example, if from any considerations it is known that the un- 
known function is close to the power, you can use the method of Sommerfeld 
(?). Its essence is to replace a segment of the power series 

f(x) = 1 + a\x + a 2 x 2 + . . . , (153) 

to the function 

f{x) « (1 + Axy. (154) 

Expanding expression (|154p in a Maclaurin series and comparing coefficients 
of this expansion with (|153p , one obtains 

a\ — 2d2 a\ 

A — ; i 1 — — o 

a\ af — 2d2 

Numerical approaches also can be used for the construction of AEFs. For 
example, in paper by ? a computational technique for matching limiting 
asymptotics is described. 

Sometimes it is possible to construct so called composite equations, 
which can be treated as "asymptotically equivalent equations". Let us em- 
phasize, that the composite equations, due to ?, can be obtained in result 
of synthesis of the limiting cases. The principal idea of the method of the 
composite equations can be formulated in the following way (?, p. 195): 

1. Identify the terms in the differential equations whose neglect in the 
straightforward approximation is responsible for the nonuniformity. 

2. Approximate those terms insofar as possible while retaining their es- 
sential character in the region of nonuniformity. 

Let's dwell on the terminology. Here we use the term "asymptotically 
equivalent function". Other terms ("reduced method of matched asymptotic 
expansions" (?), "quasifractional approximants" (?), "mimic function" (?, 
p. 181-243) also used. 

5.4 Example: Schrodinger Equation 

For the Schrodinger equation ([60]) with boundary conditions (|6T|) previ- 
ously we obtained a solution for the exponent, little different from the two 
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([68]) . In (?) the following asymptotic solutions for TV — » oo is obtained: 



(155) 



Using (|68|) and (|155)h we construct AEF 



TV 



TV + 1 



4(2iV + a)^ 



(156) 



where a = 7r 2 r(1.25) — 2 « 6.946. Numerical results are presented in Table 
5. It is evident that formula (|156p gives good results for all the values of TV. 



TV 


Eq numerical 
(Boettcher and Bender, 
(1990) 


Ea.inmi 


Error,% 


1 


1.0000 


1.0 





2 


1.0604 


0.9974 


5.9364 


4 


1.2258 


1.17446 


4.1882 


10 


1.5605 


1.5398 


1.33 


50 


1.1052 


2.1035 


0.079 


200 


2.3379 


2.3376 


0.006 


500 


2.4058 


2.4058 


^0 


1500 


2.4431 


2.4431 


^0 


3500 


2.4558 


2.45558 


^0 



Table 5. Comparison of numerical and analytical results of the energy levels 
for the SchrAudinger equation. 



5.5 Example: AEFs in the Theory of Composites 

Now let us consider an application of the method of AEFs for the cal- 
culation of the effective heat conductivity of an infinite regular array of 
perfectly conducting spheres, embedded in a matrix with unit conductivity. 
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? have obtained the following expansion for the effective conductivity (k): 

3c 



(k) 



1 



-1 



1 + a 2 c s 

~, i 

1 — ascs 



a^c 3 + a^c° 



a^c 3 



O(cf) 

(157) 

where c is the volume fracture of inclusions. Here we consider three types of 
space arrangement of spheres, namely, the simple cubic (SC), body centered 
cubic (BCC) and face centered cubic (FCC) arrays. The constants a$ for 
these arrays are given in Table 6. 





ai 


«2 


«3 


CLa 


a 5 


a 6 


SC array 


1.305 


0.231 


0.405 


0.0723 


0.153 


0.0105 


BCC array 


0.129 


-0.413 


0.764 


0.257 


0.0113 


0.00562 


FCC array 


0.0753 


0.697 


-07.41 


0.0420 


0.0231 


9.14- 10~ 7 



Table 6. The constants ai, . . . , clq in Eq. (|157|) . 



In the case of perfectly conducting large spheres (c — » c maxi where c max 
is the maximum volume fraction for a sphere) the problem can be solved 
by means of a reasonable physical assumption that the heat flux occurs 
entirely in the region where spheres are in a near contact. Thus, the effective 
conductivity is determined in the asymptotic form for the flux between two 
spheres, which is logarithmically singular in the width of a gap, justifying 
the assumption (?): 



(k) = -M 1 ln X -M 2 + 0( X - 1 ), 



(158) 



where x — 1 — ( ) 3 is the dimensionless width of a gap between the 

Y Cmax J 

neighboring spheres, \ ~~ ^ f° r c ~~ ^ c m, Mi = 0.5c max p, p is the number of 
contact points at the surface of a sphere; M 2 is a constant, depending on 
the type of space arrangement of spheres. The values of Mi, M 2 and c max 
for the three types of cubic arrays are given in Table 7. 

On the basis of limiting solutions (|157p and (|158[) we develop the AEF 
valid for all values of the volume fraction of inclusions c G [0, c max ]: 



(k) 



Pi(c) 



m+1 
P 2 C 3 



P 3 In x 



Q(c) 



(159) 



Here the functions -Pi(c), Q(c) and the constants P 2 , P3 are determined as 
follows: 
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Mi 


M 2 




SC array 


tt/2 


0.7 


tt/6 


BCC array 


VStt/2 


2.4 




FCC array 


OV^tt 


7.1 





Table 7. The constants Mi, M 2 and a t 




0.1 0.2 0.3 0.4 0.5c 

Figure 10. Effective conductivity (k) /k m of the SC array vs. volume 
fraction of inclusions c. 

<blk m 

14 
12 
10 
8 
6 
4 
2 

0.1 0.2 0.3 0.4 0.5 0.6 C 

Figure 11. Effective conductivity (k) /k m of the BCC array vs. volume 
fraction of inclusions c. 
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<hik m 




0.1 0.2 0.3 0.4 0.5 0. 6 0.7 c 



Figure 12. Effective conductivity (k) /k m of the FCC array vs. volume 
fraction of inclusions c. 



Q(c) = 1 — c — aiC3° , Pi(c)=^aics, P 2 = for n = 1, 

) + Q( 

The AEF (|159p takes into account leading terms of expansion (|157p and 
leading terms of expansion (|158[) . Coefficients are: 

1 Q(c m ax)Mi Q(c ma ^)Mi 

a = 1, «3 = ^ , «10 = Oil To 7 

OCmax 10^3 



Q(c max )M-_ 



j = 1,2, ...,m- lm, , j^3,10. 



The increment of m and n leads to the growth of the accuracy of the 
obtained solution ((159)) . Let us illustrate this dependence in the case of SC 
array. We calculated (k) for different values of m and n. In the Figure 10 
our analytical results are compared with experimental measurements from 
? (black dots). Details of these data can be found in ?. Finally, we restrict 
m = 19 and n = 2 for all types of arrays, as they provide a satisfactory 
agreement with numerical datas and a rather simple analytical form of the 
AEF ([T59D. 

Numerical results for the BCC and the FCC arrays are displayed in 
Figures 11 and 12 respectively. For BBC array the obtained AEF (|159j) is 
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compared with the experimental results from ? and ?. For FCC array the 
experimental data are not available, therefore we are comparing with the 
numerical results obtained by ? using the Rayleigh method. The agreement 
between the analytical solution (|159|) and the numerical results is quite 
satisfactory. 
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